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ABSTRACT 


A general method for steady and oscillatory/ subsonic and supersonic, poten- 
tial linearized aerodynamic flow around complex configurations is introduced. By 
applying the Green function method, a linear integral equation relating the unknown, 
small perturbation potential on the surface of the body, to the known downwash is 
obtained. The surfaces of the aircraft, wake and diaphragm (if necessary) are divided 
into small quadrilateral elements which are approximated with hyperboloiaal surfaces. 
The potential and its normal derivative are assumed to be constant within each element. 
This yields a set of linear algebraic equations. The coefficients are evaluated 
analytically. By using Gaussian elimination method, equations are solved for the 
potentials at the centroids of elements. The pressure coefficient is evaluated by the 
finite different method. The lift and moment coefficients are evaluated by numerical 
integration. Numerical results are presented. Applications to flutter are also included. 
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SECTION 1 


INTRODUCTION 

1.1 GENESIS OF THE PROBLEM 

The structural design of aircraft requires the evaluation of the aerodyncmic 
loads. Aircraft structures are flexible and this flexibility is responsible for various 
types of aeroelastic phenomena. For instance, the static aerodynamic loads depend upon 
the geometry of the surface of the aircraft, while the geometry depends upon the aero- 
dynamic loads. Therefore, it is necessary to evaluate these loads by an iterative 
process. In addition, the structural design is subject to aeroelastic constraints such as 
divergence and flutter. These constraints require accurate evaluations of steady and 
oscillatory loads. Again, the optimal design requires an iterative process, since the 
resizing of the structure changes the stability boundary. These iterative processes are 
best performed by computer analysis. This implies that considerable advantages are 
obtained if the mathematical modeling of the problem is general, flexible and efficient. 

The method used in this work is based upon the theoretical formulation developed by Morino 
(References 1 and 2). A preliminary analysis is described in Reference 3. 

It may be v/orth noting that the usual methods for the evaluation of aerodynamic 
loads are the computational lifting-surface (zero-thickness wing) theories. These methods 
are efficient and flexible but not sufficiently general for the purpose of automated design 
of complex aircraft configurations. In addition, several computational methods around 
complex aircraft configurations have already been developed. These methods are general, 
but usually quite cumbersome to use and always require human intervention to define the 
suitable types of elements to be used. In other words, it may be safely stated that none 
of the other methods currently available satisfies all the requirements of generality, flexi- 
bility and efficiency necessary for automated design. Furthermore, for oscillatory flow 
around complex configurations, only techniques based upon the doublet-lattice method 
exist in subsonic range (see References 4 and 5), while no other method is available in 
the supersonic one. 
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The method used here is applicable to steady and oscillatory subsonic and 
supersonic flows (References 6 to 10) and provides an efficient, general, and flexible 
aerodynamic tool to be used for the design of aircrafts. These advantages of the 
present method are particularly relevant in the use of automated design. 

1 .2 OUTLINE OF THE METHOD 

The present method is mainly based on the Green function method. By apply- 
ing the Green function method, the equation of the velocity potential yields a linear 
equation relating the velocity potential if, at any point, , in the flow field with 
the values of fa nd its normal derivative on the surface s, surrounding the body and 
the wake. The integral equation is obtained by imposing that the value of the potential 

at P approaches the value of <f on the surface, if \ approaches a point on the surface. 

* 

In order to solve the integral equation, the surface of the aircraft and its wake is 
divided into small quadrilateral elements. Each element is replaced by a hyperboloidal 
surface defined by the four corner points of the element. The unknown f is assumed to 
be constant within each element. Therefore, the integral equation reduces to a system 
of algebraic equations with unknown at the centroids of the elements. The 
coefficients are evaluated analytically. Once the distribution of the velocity potential 
is obtained, the pressure distribution, the lift coefficient as well as the generalized 
forces can be easily evaluated. 

In Reference 3, the surface element is replaced by its tangent plane at the 
element centroids. The improvement obtained by using hyperboloidal element are that 
the formulation becomes general, flexible, simpler to use and more accurate since the 
continuity of the surface is maintained (see References 6 to 10). 

The work presented here includes further developments with respect to the ones 
presented in Reference 8. In particular, it includes refinements and applications of the 
code SOSSA ACTS (Steady and Oscillatory Subsonic and Supersonic Aerodynamics for 
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Aerospace Complex Transportation Systems). By using this program, several numerical 
results have been obtained. These include comparison with ones obtained by lifting 
surface theory, as well as results for finite thickness wings (in both subsonic and super- 
sonic, stecdy and oscillatory flow) and results with wing-body configuration (in both 
subsonic and supersonic steady flow). Also the evaluation of generalized forces in 
oscillatory subsonic and supersonic flow is introduced. The numerical results presented 
here include the evaluation of lift and moment coefficient. Also included are a study 
of the role of diaphragms in supersonic flow, a study of convergence in both subsonic 
and supersonic, steady and oscillatory flow and a study of the effect of the truncation 
of the wake in oscillatory subsonic flow. Finally, a preliminary application to flutter 
problem is presented. 

1.3 OUTLINE OF WORK 

In Section 2 the theoretical formulation is introduced. In Section 3, numerical 
formulation as well as numerical results for steady subsonic flow are considered. In 
Section 4, numerical formulation and results for oscillatory subsonic flow are considered. 
In Section 5, numerical formulation and results for steady supersonic flow are considered. 
In Section 6, numerical formulations and results for oscillatory supersonic flow are 
considered. Preliminary flutter analysis for a rigid model with two degrees of freedom 
is considered in Section 7. Finally, conclusions are given in Section 8. Details of the 
proof of the formulation used and a list of the computer program are given in the 
Appendices. 
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SECTION 2 

THEORETICAL FORMULATION 

2.1 INTRODUCTION 

In this section, the theoretical formulation is briefly presented. Basic flow 
equations are considered in Subsection 2.2, while the boundary condition is con- 
sidered in Subsection 2.3. Green theorems for subsonic and supersonic, steady and 
oscillator/ flow, derived in References 1 and 2, are summarized in Subsection 2.4. 
The numerical formulation is considered in Subsection 2.5. The evaluation of pressure 
coefficient and generalized forces are considered in Subsections 2.6 and 2.7, 
respectively. 

2.2 BASIC FLOW EQUATION 

The aerodynamic forces are derived from the equations which relates the flow 
properties to the geometry of the lifting bodies. For many practical applications, the 
flow may be assumed to be inviscid,isentropic, and initially irrotational . In this case, 
the flow is described by the velocity potential § . In addition, the perturbation of 
the flow is assumed to be small, so that the linearized equation may be used. 

The equation of the unsteady aerodynamic potential, is 

( 2 . 1 ) 


( 2 . 2 ) 


= 


9 V 


2 

where a is the speed of sound, V is the Laplacian operator, and 


_P 

P-t 


C. _ -2 




+ v i • v 


2 - 1 



is the total time derivative. The subscript "c" reminds that should be treated as 
a constant in order to obtain the second total time derivative. 

Consider a frame of reference such that the undisturbed flow has velocity, U® , 
in the direction of the positive x-axis. By introducing the perturbation potential (f , 
such that 

$ = t <f ( 2 *3) 

(where represents perturbation velocity which is assumed to be small with respect 
to U«*> ) and neglecting higher order terms. Equation (2.1) reduces to the linearized 
equation for the potential flow 

vV - J- , o (2.4) 

Y a .* <k* 


where a- is the free-stream speed of sound, while 


^ 


(2.5) 


d-fc z»-t 

is the linearized total time derivative. For simplicity, the following nondimensional 
Prandtl-Glauert coordinates are introduced 


x - . y ■ *4 . * - ^ * %j, 


with 


P * J l-M* 

for subsonic flow and 


( 2 . 6 ) 


(2.7) 


^ = J M 3 - I = & (2.8) 

for supersonic flow. 

2.3 BOUNDARY CONDITION 

The lifting body considered here has arbitrary shape and its motion consists of 
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small vibrations of arbitrary nature. Thus, the surface of the body is represented in 
general form by 

S (%. l.-t) = O (2.9) 

Then, the boundary condition on the body is given by 

US 


- 2f * Vf -V5 - ° 


( 2 . 10 ) 


By using Equation (2.3), Equation (2.10) yields the boundary condition for unsteady 
flow. 


• Vj = ( — + IL 




) 


( 2 . 11 ) 


Furthermore, the boundary condition at infinity is given by 0, because the 
flow is assumed to be uniform at infinity. In the following, the analysis is limited 
to steady state flow. In this case. Equation (2.11) yields the steady-state boundary 
condition 

V y ' K) 5 = " H >; (2.12a) 

Consider the subsonic case first. By using the subsonic nondimensional Prandtl- 
Glauert coordinates given by Equations (2.6) and (2.7), Equation (2.12a) yields 

'dX 


(- 

<7$ 

. ■»* ♦ 

<>5 21L + 

2& 2l \ 




* 

7V } 

* 

r 

= ( _L 

.. _2. 

6 

+ ^ 


■ + 


V 

?x a x 

* Y 

^Y 

vs. 


- (J* 

$ 

** * 

'jS 


, . 



X 

<0 X 

^ Y 

^Y 

v E 

?z 

- 

+ 


- + - 

M 

: a 



1 — M 








?x p 


(2.12b) 


where N x is the component of N along x-axis (direction of the undisturbed flow). 
Therefore, 
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The second term on the right hand side is of higher order (see Subsection C.l)*. 
Therefore, Equation (2.13a) may be approximated as 


=- _l4* (2.13b) 

Next, consider the steady supersonic flow. In this case, the boundary con- 
dition is still given by Equation (2.12a). By using the supersonic nondimensional 
Prandtl-Glauert coordinates, given by Equations (2.6) and (2.8), and following the 
same procedure to obtain Equation (2.12b), one obtains 



hi ) 

M*-l vt' 

^ is the conormal derivative. 


(2.14a) 


Similarly, by neglecting the second term on the right-hand side. Equation 


• (12.14a) may be approximated as 

d* (2.14b) 


For clarity and conciseness, the boundary condition as well as the evaluation 
of the pressure coefficient for oscillatory flow are not derived here, but are presented 
in Appendix C. 


2.4 GREEN FUNCTION METHOD 

Observe that Equation (2.4) is a second order linear partial differential 


*This is true everywhere on the surface of the aircraft except at the stagnation point. 



equation with boundary conditions given by Equation (2.11) on the surface of the body 
and (f = 0 at infinity. 


The definition of the Green function of Equation (2.4) is given as 


with 


CLpo C~t 


(2.16) 


G = 0 at infinity 
where £ is the Dirac delta function defined by 

Iff / F c *•’ f - > n ’ ‘MM *)**- ,7 > 

The solution of G are derived for instance in Reference 1 . For steady 
subsonic flow, G is given by 

^ / 


(2.18) 


Lf-V R 


with 


For oscillatory subsonic flow, G is given by 


!4 


<T=' 


<HiR 


!(*.-* +T) 


(2.19) 


( 2 . 20 ) 


where ^ (t-j - t + T) is the usual Dirac- delta function and 

T-^Cr-m 


( 2 . 21 ) 


with R given by Equation (2.19) and ft ~ hj [~ 

For steady supersonic flow, G is given by 

2 “ 5 



( 2 . 22 ) 


with 


& = - _L_ 

^ sir R 

r 7/i 

R • j (X-X./- ^ C<^- %)*+ l ? ' 

For oscillatory supersonic flow, G is given by 

<T=' (5^,-^rHjrt-t.T'j) 


(2.23) 


(2.24) 


where 


T' = 



cx-x,> ± iO 


with R given by Equation (2.23) and 




(2.25) 


2.4.1 Green Theorem for Steady Subsonic Flow 

By applying the Green function method and using Equation (2.18) the 
linearized equation of the steady subsonic aerodynamic potential flow can be 
rewritten as 

4tte ^ 

TO , 1 -*)*** (2-26) 

Zw ' 

with nondimensional Prandtl-Glauert coordinates. Where E A and 2 w are the surfaces 
of the aircraft and the wake in X, Y, Z coordinate system, while N is the outward 
unit normal to , N y is the upper unit normal to 2 W and R is the vector pointing 
from any point, P*= (X*, Y*, Z*),in rhe flow field to the dummy point on surface 
'z* or , 
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(2.27) 


r - cx-x*> i +cy-Y *)) +(z-^)^ 

Moreover 

R = I r I - (cx-x*>%£Y'Y*? 2 + c z-2*Y^ 

while E is the domain function, which is defined as (see Reference 1) 

E(p,) = 1 outside 7. A 

= 1/2 on 2* 

= 0 inside 2* 

2.4.2 Green Theorem for Oscillatory Subsonic Flow 

By applying Green function method and using Equation (2.20), the 
linearized equation of the oscillatory subsonic aerodynamic potential flow yields the 
Green theorem for oscillatory subsonic flow as 



(2.30) 



(2.28) 


(2.29) 


where X* Y, Z, and 
while $ is such that 

z/r> 


5* / N, N , R, E are defined in Subsection 2.4.1 


f cx,y,z> e TAT 

$ cx,y, 2 :> €* aCT+m, ° 


(2.31) 


with 

T = Vf. ' ui /a v p r KM /^ (2-32) 

where a> is natural frequency and K is the dimensionless reduced frequency, K - • 
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2.4.3 Green Theorem for Steady Supersonic Flow 


By applying the Green function method and using Equation (2.22) / the 
linearized equation of the steady supersonic aerodynamic potential flow can be 
rewritten as 


2TTE(p,>t^>-# -U- JS. +#<*> \C~~) Jl 


where 

H * I /or Xr-X > ((fi-Y/ +G,-Zf) /z 
/<*• Xa 'X < 

II R II is the "supernorrn" of vector R and is defined by 

li R I* = I R 0 R I ^ 

with the "super-product," © , defined by 


(2.33) 


(2.34) 


(2.35) 


0. © b - &-X b x — tfy by — CL^ b z 

and the conormal derivative is given by 


(2.36) 






- /4 0 v 


(2.37) 


while R and E are defined in Subsection 2.4.1 


2.4.4 Green Theorem for Oscillatory Supersonic Flow 

By applying the Green function method and using Equation (2.24) the 
linearized equation for the oscillatory supersonic aerodynamic potential flow yields 
the Green Theorem for oscillatory supersonic flow. 
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(2.38) 


i. 

2irE cp,> f(g) -- -f£ Zfr C»iii IIRII Jl 

A 

where <p is such that 

= ^Cx,Y ? 2>e" T = $cx,y,2; e" atr '' Mx; 

* 

while the other parameters are defined in Subsections 2.4.1 and 2.4.3. 


2.5 NUMERICAL FORMULATION 

The integral equations given in Equations (2.26), (2.30), (2.33) and (2.38) 
can be solved by replacing the surface of the aircraft and wake by a finite number of 
quadrilateral hyperboloidal elements. Thus, the integral overX becomes a summation 
of a finite number of integrals over hyperboloidal surfaces which represent the surface 
of each element on . 

The general expression for the hyperboloidal surface is 

Hwpc + fP + 

(2.40) 

; -'*7*1 

where ft , p, , ft and P 3 are defined in terms of the four corner position vectors, 
p (1,1), P (1,-1), P (-1,1) and P (-1,-1). 

p c = C p Oj Of- pCl; -*> +p6-l;'> ^pC-l;-'0 
R = 4 (p O ,0 4pO/-i;-pH,0-pC-«;-)>} 

pa = ^:Cp C!> o - p Cl, -l ) +p(-l,l> - p 

_ _ _ (2.41) 

B s qKp Ct > o - P Cl, -u -pM.O + p )J 
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Assume that the values of the velocity potential and its normal derivative within 
each element are constant. Then the integrand of each surface integral become 
purely a function of the vector, p # and the geometry of the surface element on £ . 

By using Equation (2.40), the analytical expression for each element-integrc! can be 
obtained. Finally, a system of linear equation relate computed coefficients to the 
unknown, <j> k , at the centroid of element, Z k . The system is solved for <j> k by using 
the Gaussian elimination method. 

2.6 PRESSURE AND PRESSURE COEFFICIENT 

The pressure on the aircraft is evaluated from the linearized Bernoulli 

Theorem 

P ' P„ " -P. + LJ.. j (2.42a) 

< 2 r 'v > 

Therefore, the pressure coefficient is 


2 Jy. 

ul 6-t 


(2.42b) 


For steady flow, by using nondimensional Prandtl-Glauert coordinates given 
in Equation (2.6), Equation (2.42b) is reduced to 

■ (2.43) 

p 2? 

For oscillatory flow, details on the calculation of C* are given in Appendix C. 
For all the results presented here, the pressure coefficient is evaluated by finite 
different method. 
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2.7 GENERALIZED FORCES 

The evaluation of the generalized forces, used for the flutter analysis, is 
presented in this section. The generalized forces are defined as 

@h = ~ ff? n • u h d I (2.44) 

where pH is the force acting on the surface of the body, and ££ is a prescribed 
shape-function (mode). Two rigid body modes are considered in this section: they 
correspond to the lift force and the pitch moment. For the lift force, 

u h - i (2 - 45 > 

therefore, 

L 

= f/ C Pl-Pu) Jx ‘ l y (2.46) 

where fa is the z component of the unit surface normal . The subscript 11 and 

"n" refer to the lower and upper elements. For a thin wing, Q L is evaluated as 
follows (see Equation C.31) ) 



(2.47) 
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where 


AC ;> > 0 -> L 


( 2 . 48 ) 


and A is the total projected area of the wing on x-y plane, Ji is the reference length, 
Xie (■ y ) and Xf£ are the x-component of vectors pointing to the leading 
edge and the trailing edge. 

For the pitch moment, the mode is given by 

U h " <2 - » M ) I + cr-%„) ? (2 - 49) 

and corresponds to a rigid body rotation (around the axis X = ,2 = z m ) which is 

positive if the leading edge moves downward. Therefore, 

M = ff-P" + Jz 




( 2 . 50 ) 


Note that the term with p(z - z M )n^ is negligible because both z - z M and n x are 
small. For a thin wing the pitch moment coefficient is evaluated as follows (see 
Equation C.31) ) 




= t 


fit 

JJ 

z* 

<2- 

f 

'fit 
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■x 
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f» 
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SECTION 3 

STEADY SUBSONIC FLOW 

3.1 INTRODUCTION 

In Subsection 3.2, the numerical formulation for steady subsonic flow, derived 
in Reference 6, is outlined. For completeness, the proof of this formulation is given in 
Appendix A. Numerical results are presented in Subsection 3*4. 

3.2 NUMERICAL FORMULATION 

Consider the Green theorem for steady subsonic flow. According to Equation 

(2.26) 

* ei ***'•' 8 

rr (3.1) 

with boundary condition, given by Equation (2.13b) as 



By imposing that the value of the potential at P* approaches that at a point P 

on the surface X A , if P* approaches P, the value of E (P*) on the surface is found to be 

1/2 (Reference 1), and an integral equation relating the potential on the surface X. to 

A 

its normal derivative is obtained. 

Divide the surface into N small quadrilateral elements, ^ ^ , assume 
that the values of <f> and are constant within each element, and impose that 


3-1 


Equation (3.1) is satisfied at the centroid, P^, of each element 2^ . This yields N 
simultaneous linear equations, which can be expressed in the following matrix form; 




where 



is the Kronecker delta 



(3.3) 


(3.4) 


bhK C 2TT f f 


_l_ 

R 



(3.4) 


The coefficients V\/ kK represent the influence of the wake. These are obtained as follows^: 
the wake is assumed to be composed of straight vortex lines emanating from the trailing 
edge; the wake surface is divided into strips. The coefficient Wh|C = 0, if the element 
is not in contact with the trailing edge, while for the elements in contact with 
the trailing edge 


. . . cre> 

WfrK - VA/ hl< 


(3.6) 


with 




(3.7) 
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where is the strip of the wake, bounded by two streamlines emanating from the 
element . The upper (lower) sign holds for point P^, on the upper (lower) surface 
of the aircraft. 

Since no pressure difference can exist on the wake,A.<p of the wake is con- 
stant along a streamline emanating from the trailing edge, and therefore, equal to 
at the trailing edge, or approximate^ a <£> at the centroids of the elements 
adjacent to the trailing edge. Therefore, in deriving Equation (3.3), the values of 
A on the wake are approximated with the values of at the centroids of the 
elements adjacent to the trailing edge. 


3.3 ANALYTICAL EXPRESSIONS OF C/.K / and W** 

The integrals in Equations (3.4), (3.5) and (3.7) are evaluated analytically by 
approximating the surface of each element with a hyperboloidal surface described in 
Subsection 2.5. The derivations are given in Reference 6 and are outlined here in 
Appendix A. The expressions for ChK and bhu are given by 

(3-8) . 
(3.9) 

bhK r — 0 - I,H, 0 + _i J 


with 



If I 


(3.10) 
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(The subscript p reminds that the principal value of tan” from --2 to -2- 
should be considered) and 

Vn iij ^ n l , H fi + |*«i| 


+ J7, £n\\\\ |<u| * fa* 




UJ 


(3.11) 


9 f V 1/7* I 0/ . * v ✓» / 


if i r v«v 


where 


= P. + 5P, + 9P, + ?K,' ph 


(3.12) 


*. p, + n p, 


«* ' P* + f P» 


(3.13) 

(3.14) 


— „ A, * 

n ' ~ 

A, * A* 


(3.15) 


tie) 


Furthermore, the expression for Ia/^k may be obtained by taking the limit of in 

Equation (3.10) by letting J^j go to infinity. The derivation is given in Appendix 

(xe) 

A. 4. The expression of is then given as 


CrG) 


W hK * Leo- I„,C) -I Wi C-D 


(3.16) 
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with 


XmM ^ ) ' + 2 If h" 


• Pm" ( • i-x)( flt * *x) 

P-J • i«* Pj 


(3.17) 


lu/j (9)- t jj 'td'i p 


where 


-( fkjPu K x. x ’ Pj ) 4 ( ftuTfi )( ? n i • ^x) 

I P H| J * i)j * pj 


(*)= ) 


P» ! < p 4 4 f- )/ 2 p h 

P. • <p. -p-y, 

P-j - P- *•} p, 


(3.18) 


(3.19) 

(3.20) 

(3.21) 


’ In the above equations, P+ and P- are the two corner points of the element,^ , which 
are on the trailing edge. (See Figure 2b). Equation (3.3) can be solved numerically to 
yield the values of <j?^ . Once the <j?j, are obtained, the pressure can be evaluated from 
the linearized Bernoulli Theorem, given by Equation (2.43). 


3.4 NUMERICAL RESULTS 

Numerical results obtained for steady subsonic flow are presented here. In 
Subsection 3.4.1, the results are compared with those obtained by the lifting surface 
theory. Comparison with experimental results are considered in Subsection 3.4.2. 
Analysis of convergence is considered in Subsection 3.4.3. The evaluation of the lift 
coefficient is given in Subsection 3.4.4. Finally, the pressure distribution coefficient 
of a wing-body configuration is presented in Subsection 3.4.5. 
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3 . 4.1 Comparison With lifting Surface Theories 

Comparison with various lifting surface theories are presented in Figures 4-7. 
Figure 4 shows the lift distribution per unit angle of attack, C. ^ , for a rectangular wing 
with AR = 1 and M = .2. The results, obtained with NX = NY = 7 are compared with 
the ones of Cunningham 11 and Kulakowski and Haskell. 12 Figure 5 shows the distribu- 
tion of Cj( for a tapered swept wing with aspect ratio AR = 3, taper ratio TR - .5, d - 5 J 
sweep angle along quarter chord A 1/4 = 45° and Mach number M = .8. The results 
obtained with NX = NY = 7 are compared with ones of Cunningham and Kolbe and 
Boltz. 1 ^ Figure 6 shows the distribution of the section lift coefficient per unit angle of 
attack, Cj^, for a rectangular wing with aspect ratio, AR = 4 and Mach number 
M = .507. The result, obtained with NX = NY = 7 and 10, are compared with the ones 

byYates. 1 ^ 

All these results are in good agreement with the existing one. Thickness ratio 
used in evaluating the above results is chosen to be 0.1%. The reasons to choose this 
value is that, as shown in Reference 3, the solution withl> 0.1% is a good approxi- 
mation for the zero thickness solution and that values of thickness ratio smaller than 
0.01% may cause elimination of significant figures. 

3.4.2 Comparison With Experimental Results 

Figues 7, 8a and 8b present a comparison with the experimental results of 
Lessing, Troutman and Menees. 15 The results are relative to a rectangular wing with 
aspect ratio AR = 3. The airfoil consists of a biconvex circular arc section with 5 
percent thickness and with sharp leading and trailing edges. Figure 7 shows the pressure 
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distribution on the upper and the lower surfaces, £^ u an d , versus x/c at 2y/b = 

0., 0.5, 0.7 and 0.9 for the above mentioned wing witho(= 0° and M = 0.24. 

. Figure 8a shows the lift distribution,^ , versus x/c at 2y/b = 0, 0.5, 0.7 and 0.9 for 
the above mentioned wing with c<= 5° and M = 0.24, while Figure 8b shows the 
pressure distribution on the upper and the lower surface, C^ u and versus x/c at 
2y/b = 0, 0.5, 0.7 and 0.9 for the above mentioned wing with = 5° and M = 0.24. 

It should be mentioned that the results for Figure 8a are obtained by taking the 
advantage of antisymmetry. If the velocity potential,^, is separated into tv/o parts, 
symmetric part,^ , and antisymmetric part, , then the difference of the velocity 
potential on the upper and lower surfaces of the wing is twice the value of ^ while the 
velocity potential on the upper surface of the wing is t* and the velocity 

potential on the lower surface is 4? = ' fes • Tbe VQ l ue of can be obtained 

by taking advantage of antisymmetry while the value of can be obtained by taking 
advantage of symmetry. This is mainly because the doublet, source and wake 
integrals are nearly independent of angle of attack, eA . Therefore C ^ , bf,x and 
Whk can be evaluated at —0, while the onlyc4 -dependent term in Equation (3.3) 
is the linearized boundary term, , which can be separated into symmetric and 




antisymmetric as follows: Set 


= - n x - -| W cos*< ±('r\ z ) srn*0 

"3 |SJ s e<-o 

- - (nx)j so t 
and as mentioned before, 

+ = i +« 


& 


(3.22) 


(3.23) 
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where (h*) andC^) are the x-component and z-component of the surface unit normal . 
The upper (or lower) sign holds for the upper (or lower) surface. A similar example of 
the symmetric and antisymmetric problem is given in Section B.5. It may be worth 
noting that the results for Figure 8b are obtained by directly solving <Jj, or with C/,^ , 
and evaluated at ck = 5°. However, the lift coefficient obtained by 
subtracting Q ^ from is consistent with results obtained for Figure 8a. 

3.4.3 Analysts of Convergence 

Figure 9 presents a study of convergence of the values of in terms of the 
number of elements. The wing is the same as that of Figure 7 with <k — 0 and M = 0.24. 
In the results shown in Figure 9, there is a sharp change of the values near the leading 
edge and the wing tip, therefore, in order to increase the rate of convergence, 
smaller elements along the leading edge and the wing tip are used. For this purpose, 
the coordinates of the nodes (x, y) are defined by 

Xr - Cihr)* ' - °> 1 > 2 — ^ 

(3.24) 

-=r ^ 

j = 0 , 1, '2 N Y 

Result shown is the distribution of along y = 0. It may be noted that 100 elements 
on the whole wing (i . e. , NX = NY = 5) are sufficient for convergence. 

3.4.4 Total Lift Coefficient 

The total lift coefficient for steady flow is evaluated from Equation (2.47) 
without the imaginary terms, i.e.. 
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(3.25) 




where <jy and ^ are the extrapolated values of the velocity potential along the lower 
and upper trailing edges. Figure 10 shows the values of Cl of a triangular wing with 
different aspect ratios in incompressible flow. The result obtained at "C = 0.1% are 


compared with the analytical ones. 


16 


3.4.5 Wing-Body Configurations 

In Figure 11, results for the pressure distribution coefficient of a wing-body 

configuration in steady subsonic flow are compared, with the results presented by 

1 7 

Labrujere, Loeve and Slooff. The results were obtained for M = 0 and a rectangular 
midpositioned wing with chord c = 1, span b = 6, thickness ratio ~C= 0.9, and angle 
of attack cAj = 6°. The body is at zero angle of attack and is composed of a forebody 
with length — 2 and radius 


f = 0.5 - 0,125 


(3.26) 


a midsection of length — 1 and radius r — 0,5 and an aftbody of length — 9, ana 
radius r = 0.5. It is worth note that the extra long aftbody is used to simulate the wcke 
emanating from the circular fuselage. The number of elements is 200 on the whole 
configuration (NX = 5, NY = 5 on the wing, NX = 2, NY = 3 on the forebody, NX “ 
NY = 3 on the midsection, and NX = NY = 3 on the aftbody). Figure lib presents 
the distribution of cb - along three circumferential stations (see Figure lib) for the 
same wing-body configuration. 
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SECTION 4 

OSCILLATORY SUBSONIC FLOW 


4.1 INTRODUCTION 

In subsection 4.2, the numerical formulations for oscillatory subsonic flow 
are considered. The boundary condition and the pressure coefficient are considered 
in Appendix C. In subsection 4.3, numerical results obtained are presented. 


4.2 NUMERICAL FORMULATION 

Consider the Green Theorem for oscillatory subsonic flow, given by Eq. (2.30) 




WAR 


,WAR 






a - . 




s„ • f? 

with boundary condition given by Equation (C.24) 


(4.1) 


2-L =k[ ( ‘i v'V + 1 ^2 \ _ 




(4.2) 


«t4 ? <?X 

( — ' 

where 2 . (X,Y) is the vibration mode defined by Equation (C.18) 

By imposing that the value of the velocity potential at*P* approaches that 
ot a point, P, on the surface if P* approaches P, the value of E(K) is found to 
be 1/2 (Reference 1), and an integral equation relating the potential on the surface 
to its normal derivative is obtained. Noting that no pressure difference can 
exist on the wake, and using Equation (C.32) 


one can conclude that 0^)6^ is constant along a streamline emanating from 
the trailing edge and equal to the value at the trailing edge 

, 4 . 4 , 

In Equation 4.4, the value of ^6* ^/approximated with this value at 

the centroid of the element^ adjacent to the trailing edge. In Equation 4.4, X « 
represents the X component of the position vector of the centroid of elements, 

, adjacent to the trailing edge where the wake starts to develop. Using Equation 
(4.4) and applying the same procedures used for the steady flow, one obtains 




,A) 


(4.5) 


where 




r**r» 


(4.6) 


- 1*8 




%.< V' 

and for elements not adjacent to the trailing edge, W hK =0, otherwise. 


(4.8) 

2k <* lh 

where X K is defined previously. £ ^ is the strip of the wake, bounded by two 
streamlines emanating from the element J The upper (or lower) sign holds for 
point, p^, on the upper (or lower) surface of the aircraft. Using the following 

equation 


(4.7) 


~i A K 


vAi - (l 4 & 


( 4 . 9 ) 
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the approximate evaluations of C hK and b hk are obtained by replacing the exponen- 
tial term with its value at the centroid of the element. This yields 


A -t&K c 

C-ly\Z ~ C ^ 1 ^ Oi K 


( 4 . 10 ) 


b h y = e b hK (4ji) 

where C hK and b h(< are given by Equations (3.8) and (3.9) whereas R c is the distance 
between the centroids of elements^ and 1 ^ For the wake contribution one obtains 
, if the element ^ is not in contact with the trailing edge, while, for the 
element 5, K in contact with the trailing edge, dividing the stripX^ into small elements, 

one obtains 

A ACT*) NW + 

J p o + .^-fe)a j (4 .i2) 

iP‘ 

where, j refers to the jth wake element emanating from the element , C^. is given 
by Equation (3.8) R. Is the distance between the centroids of element, £ h , and the 
ith wake element; NW is the number of elements of the truncated wake emanating 

from (A more refined analysis is obtained by taking the limit of the present 

one when NW tends to infinity) . 


4.3 NUMERICAL RESULT 

Results for oscillatory subsonic flow are presented here. In Subsection 4.3.1 the 
results are compared with the experimental ones. In Subsection 4.3.2 the analysis of 
convergence is considered. In oscillatory flow, instead of using semi-infinite wake, 
trancated wake is used. Therefore a study of the wake length LW as well as the 
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analysis of the convergence due to the number of wake elements are considered 
in Subsection 4.3.3. 

4.3.1. Comparison with Experimental Results 

Figure 12 presents a comparison with the experimental results of 
15 

Lessing, Troutman and Menees . The results are relating to a rectangular wing 
with aspect ratio AR = 3. The airfoil consists of a biconvex circular-arc section, 
5% thickness with sharp leading and trailing edges. The wing considered has the 
same geometry as the wing considered for Figures 7 and 8 and is oscillating in 
bending mode described by 

2 = ,i?04-3 |-r^l *1* 1.70255 l~jfl ~ | -Xl 4-.252S7 |t^ j 

b (4.13) 

with 

M (M * 0.24- 

The values showed are the real and imaginary part of the distribution 
ofC L along *^ 3 0and *^*0.5 • They are obtained with NX=NY=7, NW=20, and 
LW-2 C. The node coordinates are given by Equation (3.22). 


4.3.2 Analysis of Convergence 

A convergence study of the problem considered in above subsection 
are presented in Figures 13a and 13b. The values showed are the real and 
imaginary parts of the distribution of along and 2y/b=.1328 
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respectively. They are obtained with different numbers of wing elements, 

NX=NY= 5, 6,, 7, while the number of wake element emenating from the 
trailing edge and the wake length, LW, are fixed such that NW=30, LW=3.5C. 
Using 64 elements on the whole wing (i.e. NX=NY=4) is sufficient for convergence. 


4.3.3 Analysis of Wake Effect 

The wake effect of the oscillatory subsonic flow is evaluated from 
A 

Equation ( 4 . 12 ). The values of on the number of wake elements and 

the wake length. Note that, C . in Equation (4.12) is a doublet integral which 

hi 

represents the solid angle. The value of the solid angle decreases as the distance 
between the control point and wake element increases. Therefore, the wake effect 
is dominated by the wake elements closest to the trailing edge. For this reason, 
it is necessary to determine the effective wake length. Also after the effective 
wake length is determined, it is necessary to study the convergence due to the 
number of wake elements. Figure 14 shows the distribution of C ^ , along 
7^=0.13265, of the same problem considered for Figure 12a . The results 
are obtained with different number of NW-5 , 10, 20, 30, 40, while using 196 
elements on the whole wing, i .e., NX=NY=7, and using four times chord length 
for the wake length, or LW -4.0C. The length of the wake elements varies from 
0.8C to 0.1C from the figure shown. At least NW=20 should be used for converg- 
ence. Figure 15 shows the distribution of 'c^r along 2)/^ =0.13265, of the same 
problem considered for Figure 12a. The results are obtained with different lengths 
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of wake LW — 0.5C, C, 2C,3C and 4C. The length of the wake element is 

the same for each run, or (aXJ w =LW/NW = 0. 1C. The numbers of NW for each 
run are 5, 10, 20 and 40 respectively. From the figure shown, it can be con- 
cluded that at least, LW=2C, should be used as effective truncated wake length. 


4.3.4 Results of Generalized Forces 

Figures 16a and 16b shows the absolute values and the phase angles of 
Cl and for a delta wing in incompressible oscillatory flow , as functions of the 
reduced frequency, K= The wing with AR=4 and -f = . 5% is oscillating in 

pitch about the axis X=X H . The vibration mode is described by with 

X., =0.5 (X + X ). The results are compared with the experimental, as well 

M LE TE 

18 

as, the numerical results by Laschka 

Figures 17a and 17b shows the absolute values and the phase angles of 

C, an< ^ C as function of the Mach number which varies from 0.0 to 2.5. (Section 
L M 

5 presents the formulation for the supersonic range) for a rectangular wing with 

AR=2 and "C ^0.001 . The vibration mode is the same as that of the above problem, 

and the reduced frequency is set at K=1 . The results are compared with the ones 

18 

obtained by Laschka . Figures 18a and 18b show the thickness effect on the 
problem of Figure 17. The results are obtained with two thickness ratios 
=0.001 and TT=0.05. The Mach number varies from 0.0 to 2.5. 

Figures 19a and 19b show the absolute values and the phase angles 

. 

of C, and C. . as functions of the reduced frequencies for a rectangular wing 

L M 
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with AR=2 and ~C =0.001 . The wing in incompressible flow is oscillating in plunge. 

The vibration mode is described by e ,u,t . The results are compered with the ones 

by Lawrance and Gerber^ as well as Laschka . 

Figures 20a and 20b shows the absolute values and the phase angles of 

and^ of the same problem as that given for Figures 19a and 19b, while the 

L M 

reduced frequency is set at K=1 .0 and the Mach number varys from 0.0 to 2.5. 

18 

Results are compared with the ones of Laschka . 




SECTION 5 


STEADY SUPERSONIC FLOW 

5.1 INTRODUCTION 

In Subsection 5.2 the numerical formulations for steady supersonic flow, 
derived in Reference. 7, is outlined. For completeness, the derivation of these formula- 
tions and more detail formulations are summarized in Appendix B. In Subsection 5.3 
wings with supersonic leading edges and diaphragms for wings with subsonic leading 
edges are discussed. In Subsection 5.4 numerical results are presented. 

5.2 NUMERICAL FORMULATIONS 

Consider the Green theorem for steady supersonic flow. According to Equation 

(2.33) 


27TE = „ 4- J £+ tii-2- (5.1) 

1 1 u IIRII A {£' ^ 


with the boundary condition given by Equation (2.14a) 


*£ = _ k\ (± + Jil 




= -M^-L + JTL ®3L) 

x £ M J - 1 'i'i. ' 


(5.2) 


Applying the same procedure for subsonic flow, one obtains the following 
equations, expressed in matrix form, 


( Hk) fie ‘ ” C khj 


(2±) 


(5.3) 
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where 


11^,1 


and 


(5.4) 


bhK "^ *{[ ifr d ^-* 


P**Ph 


(5.5) 


For hyperboloidal elements completely Inside the Mach forecone, anc ® 

t>#, (c are still given by 


ChK = Ip c I , I ) - J v c\ r I) -X/-I, I ) + L c-b-i ) 


(5.6) 


b bK = U Cl . <> - I 5 CI, -O-IsC-b I) +l 5 (~b-l) 


(5.7) 


For elements completely inside the Mach forecone Ip and 1^ are given by 


and 


i.c 

s r ''f 1 h %<«* 




(5.8) 


(5.9) 
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with 


h = -C jp.tf,-,) 


(5.10) 




(5.11) 


M.l) = 

8 ” 1 iif ii }•»,/«; / 


(5.12) 


where q, a^, 02 / and n are given by Equations (3.12) through (3.15) while 

1 1 W 0 


11 til 


f® *j 


a j e °j =0 


(5.13) 


- -~Sr~ Sin’' f-t^L) J.o7.<<, 

IKjll 1 >i|* Ajii 2 i J 

For elements partially inside the Mach forecone Ip and T^-j, I^/ I 53 are still 
given by Equations (5.8), (5.10), (5.11) and (5.12) if the comer point is inside the 
Mach forecone. However, if the corner point is outside the Mach forecone, integrals 
in Equations (5.4) and (5.5) involve singularity. Details of these kinds of integrals 
are discussed in Subsection B.4. For elements completely outside the Mach forecone. 


H = 0 


(5.14a) 
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Therefore 


“ thK ' ° 


(5.14b) 


Note that the wake effect has not been considered for the results presented 
here, because the geometry of the lifting body under consideration includes only wing, 
nose, or middle-section fuselage. Therefore, for each element, the wake is always outside 
the Mach forecone, i.e., WhK = 0. 

In addition, by neglecting the higher order terms, the boundary condition is 
given by (see Equation 2.14b) 


£ ■ - "■/(> 


5.3 DIAPHRAGMS 

There are three categories of wing geometry in supersonic flow: 1) supersonic 
leading edges (the elements on the upper or the lower sides are influenced only by 
the elements on the same sides, therefore, the integration in Equations (5.4) and (5.5) 
are around the upper or lower surface only for the upper or lower elements); 2) mixed 
supersonic-subsonic leading edge (for this kind of wing, diaphragms have to be used to 
separate upper and lower sides of the aircraft; then, the problem can be handled as the 
case of supersonic leading edge); 3) subsonic leading edge (for this kind of wing, 
results obtained with or without diaphragms are the same [see Subsection 5.4.5]; 
however, for small thickness ratio, diaphragms are still needed to avoid the elimination 
of significant figures). 

If diaphragms are used, both values of the velocity potential and its normal 
derivative, and ^ , of the diaphragm element are unknown. However, two 
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equations are obtained for each element of the diaphragm: one relates <p p and 

of the diaphragm element, , to the upper geometry of the aircraft and the 
upper side diaphragm, the other one relates the same quantifies to the lower geometry. 
Therefore, the total number of equations is equal to the total number of unknowns. 
Details of the numerical procedures are given in Subsection B.5. As mentioned before, 
for wing with subsonic leading edge, procedure with or without diaphragm can be used. 
If the procedure without diaphragm is employed, then the integral equations are solved 
the same way as for the subsonic case (solutions are obtained by solving Equation (5.3) 
However, if the procedure with diaphragm is employed, solutions are obtained by 
solving Equations (B.41) or (B.46). 


5.4 NUMERICAL RESULTS 

Results for steady supersonic flow are presented here. In Subsection 5.4.1, 
results are compared with the experimental ones. In Subsection 5.4.2, analysis of 
convergence isstudied. Problems ofairfoilsin conical flow are studied in Subsection 
5.4.3. The role of diaphragms has already been discussed in Subsection 5.3. In order 
to show that diaphragms are not necessary for geometries with subsonic leading edge, 
two results obtained with and without diaphragms for a circular cone in conical flow 
are compared in Subsection 5.4.4. Finally, in order to study body interference and 
also to show the generality of the present method, result of wing-body configuration is 
considered in Subsection 5.4.5. 

5.4.1 Comparison With Experimental Result . 

Figure 21 shows the distribution of the pressure coefficient Cp on the lower and 
upper surfaces for a rectangular wing with aspect ratio AR = 3. The airfoil consists of a 
biconvex circular-arc section, 5% thick, with sharp leading and trailing edges. The 
results are obtained forc^ = 0* M = 1 .3 and NX = NY = 7. Figure 22a shows the 
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distribution of the lift coefficient for the same wing considered in Figure 21 . The 

results are obtained for = 5°, M = 1 .3 and NX = NY = 7. These results are com- 

15 

pared with the ones obtained by Lessing, et al. In Figure 22b, the pressure 
distributions on the upper and the lower sides of the wing are shown. It should be 
mentioned that results in Figure 22a are obtained by taking the advantage of anti- 
symmetry as in the case described in Subsection 3.4.2. 

5.4.2 Analysis of Convergence 

Convergence study of the problem considered in Figure 21 is presented in 
Figure 23 which shows the distribution of the velocity potential along 2y/b = 0 for 
different numbers of elements. The curves are obtained with NX = NY = 5, 6 and 

7. It is shown that 144 elements on the whole wing, i.e., NX = NY = 6, are 
sufficient for convergence. 


5.4.3 Delta Wings 

Figure 24 shows the distribution of life coefficient per unit angle of attack for 
a delta wing with supersonic leading edge and 

% -A- r 1,2 

where _A- is the sweep angle of the leading edge. The results obtained with NX = 8, 

NY = 12 and M = 1 .2 are compared with the exact solution which is given by (Reference 

20 ) 


A 


Cp c — M_ 



f-VnV 
VM -Y 


+ Cos"' 


jjf Vwif "I 

w+ar s 


(5.16) 


where 
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Figure 25 shows the distribution of the li ft coefficient per unit angle of 
attack for a delta wing with subsonic leading edges and m = 0.8333. The results 
obtained with NX = NY = 7 and M = a[~ 2 are compared with the exact solution, which 
is given by (Reference 20) 


A C - 4- ft V _J (5.17) 

^ 

where K = r\ 1 - and E is the elliptic integral of the second kind. It may be worth 

noting that there are two types of elements — grids which can be used for delta wings 

(see Figure 2c). For supersonic delta wings, the flow is conical, therefore, the 

element-grid shown in Figure 2c (b) appears to be the more natural one. In addition, 

the use of this grid implies the use of a completely general quadrilateral element. Let 

and tj be the generalized coordinates such that ^ has same direction as ^ while 

the direction of ^ passes the origin and pointed downward. The derivative of <£> with 

respect to X is obtained as follows: consider the directional derivatives and — 

dr) 


J<{> 

d 7 


tiL- . -it-T \ = 

l * 4 J ' ||| 




l?l 


where 


dtp 


v* J Y* 3YJ xVy> 


•,(5.18) 

(5.19) 


^ x * y j 

n - Y J- 


(5.20a) 

(5.20b) 
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Therefore 


?k _ /At \ 

' ^ dj ) 




Jr+r ( A ) 
x l V 


jl 

X 


(5.21) 


The values of (^t-) and ( can be evaluated by the finite difference method. 
K o-j ' 'dr) * 

By using Equation (5.21), the values of C ^ are obtained. The results obtained above 
by applying this scheme are very accurate. The results obtained for Figure 24 are 
especially accurate (up to the fifth significant figure). 


5.4.4 Circular Cone (Analysis of Diaphragm Existence) 

For most cases of supersonic flow considered here, diaphragms seem to be 

necessary to avoid determinant singularity. However, in this Subsection, an example 

is given to show that diaphragms are unnecessary, for example, for bodies with subsonic 

leading edge. This example is related to a circular cone with base radius r = 1, length 

I = 5, and M = 2.0. One is obtained with diaphragm by solving Equation (B.46), while 

the other one is obtained without diaphragm by solving Equation (5.3). The results are 

in good agreement. The pressure coefficient of the body is evaluated from the slope of 

curve. This value is evaluated as 0.112 which is consistent with the exact solution 

given in Reference 21 (Figure 17, p. 187, C D *= 0. 12 for M = 2, 0 = tan“^ ( — ) - 11°). 

p c 5 


5.4.5 Wing-Body Configuration 

The present method is sufficiently general to be applie. for any complex con- 
figuration. Presented here is an example of this application. A wing-body combination 
in supersonic flow with M = 1 .48 is considered in Figures 27a, 27b, and 27c . The 
geometry of the configuration is shown in Figure 27c. 

The combination is composed of a wing with chord c= 3, span 5=9, thickness 
XT = 5%, a forebody of length = 6.0, and radius varying from 0 to 0.75 linearly, and 
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a midsection of length L M 3 *° and radius r — 0.75. Wake and aftbody are not 

considered. The angle of attack of the wing is — i .92, while the angle of attack 

of the body is zero, = 0. The results are obtained with 580 elements on the whole 

configuration (NX = NY = 10 on the wing, NX = 5, NY = 3 on the nose, NX = 10, 

NY = 3 on the middle section). In Figure 27a the distribution of on the 

22 

wing section is presented, and the results are compared with the ones by Nielsen 

23 

and Woodward, Tinoco and Larsen, while in Figure 27b the distribution of C ^ /j w 

on the fuselage is shown. 
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SECTION 6 


OSCILLATORY SUPERSONIC FLOW 
6.1 INTRODUCTION 

In Subsection 6.2, the numerical formulation for oscillatory supersonic flow is 
considered. The boundary conditions and the pressure coefficient are considered in 
Appendix C. The role of the diaphragm in oscillatory supersonic flow is the same as 
in steady supersonic flow. It has been already discussed in Subsection 5.3, therefore 
if is unnecessary to repeat it here. In Subsection 6.3, numerical results obtained are 
presented. 


6.2 NUMERICAL FORMULATIONS 

Consider the Green Theorem for oscillatory supersonic flow. According to 
Equation (2.38) 


21TE (p*;4>(p*)= — C .06 (JL 

' 7 ‘ 7 HR II 


(6.1) 


+ -£ k (7T. 


HR l| 


with the boundary condition, given by Equation (C.28) 


21 


hl z (r k ST + 




\ i-AMX 

) e 


( 6 . 2 ) 


where Z (X,Y) is the vibration mode defined by Equation (C.18). Applying the same 
procedure used for subsonic flow, one obtains the following equations, expressed in 
matrix form 
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(6.3) 


C£j 


r* 




where 


C hK'~(-v ff v^im 005 ^ 1 |R| 0)^J 


(6.4) 


fie a Ph 


and 



(6.5) 


The wake effect for supersonic flow is not considered here since only 

A 

supersonic-traiiing-edge wings are presented here. Therefore W},k = 

Noting that 

- - a [o*> {/i upii) + ^ ii^n (-°- (6,6) 

- J^Oo-o (_Ti- HR ll) -f' l|ff II <aiv) {sx. |IR||)) ~~^c. ( ]j|j) 

A A 

the approximating evaluations of Q ancl b hK are obtained by replacing the sine and 
cosine terms with their values at the centroid of the element. This yields 


A f 


CfcK = [coo(a |||? c l|)+uQ||Rcll stn{n \\M)}C h t 


(6.7) 
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and 


bhK r COO (-^ USd 1 ) bhK 

where anc ^ bf^ are 9' ven by Equations (5.6) and (5,7) with Ip, 1^ given by Equa- 
tions (5.8) through (5.12) for elements completely inside the Mach forecone. For 
element partially inside the Mach forecone. Ip and I $ are still given by Equations 
(5.8) and (5.9) if the corner point is inside the Mach forecone, otherwise, they are 
given by Equations (B.26) through (B.28) for I^, Equations (B.30) through (B.32) for 
I 0 , and Equations (B.35) through (B.40) for Lg and Ip. For element completely 
outside the Mach forecone Cmc = khl< = 


6.3 NUMERICAL RESULTS 

Results for oscillatory supersonic flow are presented here. 

In Subsection 6.3.1 a wing oscillating in bending mode is considered. In 
Subsection 6.3.2, the convergence problem is considered. In Subsection 6.3.3, the 
evaluation of generalized forces is considered. 


6.3.1 Comparison With Experimental Results 
— ■ 

Figure 28 shows the distribution of absolute value and phase angle of C \ for 
a rectangular wing oscillating in bending mode 


2 * .im3 \^\ + \ •10255 1^1*- I . I ? 588 1 ^ J 3 +. 2 53% f ]-^ I (6. 9) 


with |< ^ U~ =0-1 

The airfoil, with aspect ratio AR = 3, consists of a biconvex circular-arc 
section, 5% thick with sharp leading and trailing edges. Results obtained for M = 1 ,3 
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and NX = NY = 7 are compared with the ones by Lessing, et al. It should be men- 
tioned that the sign of phase angle of obtained by the present method is opposite 
to that of Ref. 15. The author believes that this is due to the different choice of the 
sign convention. 

6.3.2 Analysis of Convergence 

A convergence study of the above problem is presented here. In Figure 29, 
the distributions of the imaginary and real parts of the velocity potential along 2y/b = 
0.5 for the problem given for Figure 28 are shown. The curves are obtained with 
NX = NY = 5, 6, 7 and 8. It is shown that 100 elements (i.e., NX = NY = 5) are 

sufficient for convergence. 

6.3.3 Results of Generalized Forces 

The typical results for the evaluation of generalized forces are presented in 

Figures 17a, 17b, 18a, 18b, 20a and 20b. In Figures 17a and 17b, a rectangular 

wing with oscillation in pitch is considered. In Figures 18a and 18b the thickness 

effect of the problem given for Figures 17a and 17b is considered. In Figure 20a and 

20b a rectangular wing with oscillation in plunge is considered. These results are 

obtained for Mach number varying from 0.0 to 2.0 and K = —7; — = 1 .0 and compared 

18 U “ 

with the ones obtained by Laschka. Details of the problems have already been given 
in Subsection 4.3.4, therefore, they are not repeated here. 
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SECTION 7 


FLUTTER ANALYSIS 


7.1 INTRODUCTION 

The previous sections present an efficient aerodynamic tool to analyze 
flutter phenomenon. In this section, a preliminary flutter analysis is presented. In 
Subsection 7.2, the mathematical model is presented. In Subsection 7.3, the Lagrange's 
equation forces are discussed. In Subsection 7.5, the evaluation of aerodynamic forces 
is outlined. In Subsection 7.6, the flutter analysis is described. In Subsection 7.7, 
numerical results are presented. 

7.2 MATHEMATICAL MODEL 

As usual, the mathematical model is established by describing the aerodynamic 
loads due to simple harmonic motion and then studying the stability of infinitesimal 
vibration of this kind of motion. The model is based on the fact that the motion at 
the flutter boundary is simple harmonic and that the flutter boundary is by definition 
the one relative to an infinitesimal oscillation around the next configuration. 
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The most typical aircraft flutter results from a coupling between 
the bending and torsion modes of a wing. A typical section (airfoil) analysis is 
considered for simplicity. This kind of flutter is considered in this section. The 
airfoil under consideration is permitted freedom to execute small vertical and 
angular displacements/ h and ck , and is constrained by two springs representing 
bending and torsion respectively. Thus, the displacement in 2 direction is 


(7.1) 

where X denotes the position of elastic axis and h is the vertical displacement of 
EA 

the elastic axis (Figure 3). 

Mathematically, it is easier to describe displacement with exponential 
time dependent e'^*’ (id real) and assume that all dependent variables are proport- 
ional to e' W *. Therefore, Equation (7.1) can be rewritten as 


W4e i, % (x -W?e i, ‘ 


w e 


(7.2) 


with 

Qj « 'h + £(X - X EA ) (7.3) 



d 




The actual quantities of W/ h and are 


(7.4) 

(7.5) 


found by taking the real part of the complex 


counterparts. 
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7.3 


LAGRANGE’S EQUATION OF MOTION 

Equation (7.1) describes the displacement of airfoil having two 
degree of freedom. The variables h and<X are two generalized coordinates. The 
motion of the unrestrained airfoil is defined by the Lagrange's Equation which 
states that 


J / *T \ '<7 » , y «-J _ 


±7 4 2 ^ 


(7.6) 


where i$ the k-th generalized coordinates, T is the kinetic energy of the 
wing segment, U is the internal strain energy of the wing segment, Q k is the k-th 
generalized force corresponding to the k-th generalized coordinate. The kinetic 
energy of the wing segment is 
*r - — ~ [W Jw 

2 s c hord 

c F + ( X *“ 

S chord 

* - (7.7) 

where Jl is a reference length, while r^ and X* are the dimensionless radius of 
gyration and static balance about the elastic axis. 


Yd. ' T(X-Kea) dm j ^ J(X Xu) di 


(7.8) 


If the stiffness of bending and torsion spring are represented by and 
respectively, the internal strain energy is given by 
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or 


Ue * m «)h \n + t >" CrJ «; ^ 


(7.11) 


where tl/^ andU/^ are the natural frequencies without mass coupling, or 




Kk 


c^, 


_ rur 


(7.12) 


(7.13) 


Combining Equations (7.7) and (7.11) with the Lagrange Equation, Equation (7.6), 
with cjj- h, q 2 =^/ Qj=Q fc ancl } f°" ow ' n 9 two equations of motion are 

obtained 


Mb +■*! + W h = $h 


(7.14) 


h + (7.15) 

Combining Equations (7.4) and (7.5) with Equations (7.14) and (7,15) yields 

^ ^ m Wh h - &h 

(7.16) 

- in w*&x<t'h - w i srxyi + wi y; (7.17) 

As mentioned above, the symbol tilde above each quantity represents the complex time* 
independent part. 


* .i 


- hi W*h " •+ ™ - &h 




7.4. 


GENERALIZED FORCES 

The generalized forces are found from the virtual work which is given 


by 


* ff* c ? > Ax * i 
~ 2 /, fh + 


(7.18) 
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where 


3 * * L 


(7.19) 


and 


^r- ff (7.20) 
Combining Equations (7.19) and (7.20) with Equations (7.14) and ( 7.15) yields 


£ hi i«W J) h ~ w K )(* ~ L 


(7.21) 


- hx + (“ ww; * J? a fZ- + W at; i Yl )° 4 s ^ ^ 


7.5 AERODYNAMIC FORCES 

r** rv 

Following is the procedure to evaluate L and M. Because of the 
assumption of small vibration, the doublet integral and source integral (see Equations 
(4.4), (4.5), (6.4), and (6.5) ) are not influenced by the displacement, W, while 
the following boundary condition is a function of W 


<b§> i , i 0 \A /\ 

^fsJ = ^ ^ X ) 

Therefore, it is legitimate to write 

C f ‘ = i w '> 


-- i(X e!^ * ^ 

* X X.Ce’ wt ) * 2m-**) &"*) 


(7.23) 


(7.24) 
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where "jL denotes the linear operator. 

Therefore, by applying Equation (7.24) to Equation (7.19) and (7.20) 
the lift force and pitch moment are found. 

Z * 1/ ^ ^ 

j'Ot- 




(7.25) 


* e ~ (-f- Clk+Z cj) 


2 
i a 


M r ~ f/Cp CX-X 6 ^t)x 

+ ? /// x - X e< ) g'' 1 " 4 ; (X -X^) Jx Jy j 

•< O 


<° £ z [ r + 

o V / ^Mh 


(7.26) 


In Equations (7.25) and (7,26) and C ^ are lift force and pitch moment 
coefficients about elastic axis of the airfoil in the motion with oscillation in plunge. 

h 

The displacement for plunge motion with -j- = f is described by 

w . I ^ 

The boundary condition is then given by 


(7.27) 


= nJ f 
'ilJ N ” 1 


• A |wt N 

ul^ €. ) 


(7.23) 


or 


= fv/ z A- £ 
<2 N 


(7.29) 



Similarly, Cu anc ^ are l'f* force and pitch moment about elastic axis of the 

airfoil in the motion with oscillation in pitch about the elastic axis with unit angle 
ampitudco The displacement is described by 


iwt 


W c OC £ 

the boundary condition is then given by 


(7.30) 


-- lJ z Cl + : " Ctt-XaO 

t 




i.e. 


^2 + «k (xdhtD 

£ 


(7.31) 

(7.32) 
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FLUTTER ANALYSIS 

Combining Equations (7.25) and (7.26) with Equations (7.21) and (7.22) 


yields 

+ e ^- c G* 


r-' 

A 




ct- 0 (7.33) 


and 


f 




_ f d* c 




By introducing 


and 



k 


(7.35) 

(7.36) 
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(7.37) 

(7.38) 


Equation: (7.33) and (7.34) may be rewritten in the following matrix form 


1 Qa Cu 


' n ’ irK*/* 

*2. 


X<*^ TiKyM 


s. 

[VI 


*• 


ry 

/ 1 

d 


- 0 


(7.39) 


side matrix must be zero. The value of A can, therefore, be evaluated. By plotting 
the values of the imaginary part of A versus different reduced frequencies, the 
flutter is found at a specific reduced frequency where the imaginary part of A is 
zero. 


7.7 NUMERICAL RESULTS 

In order to compare the numerical results with the analytical ones, a 
two-dimensional thin airfoil in incompressible, oscillatory flow was studied. For 
this airfoil, the analytical expressions of Q h / / C Mfl and ai 

given in several standard text books (see e.g. Reference 24) as 


C Lh -titiec* j? * -i) (7 ‘ 40) 

^.4,) 

^ = (7 - 42) 
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ui^+(2h+i)^(^o0yji<-iljl (7 43) 

where 

Ck = (-J, (h)+:y,w)/i (-(J'M+y.m) f ; (y,(k) -xadj (7.44) 

where J ^ and and J-Bessel functions of the first and zero orders and Y| and Y 

are Y-Bessel functions of the first and zero orders and reduced frequency, K, is 
defined as 

i< = a** 

(7.45) 

It should be noted that the purpose for this two dimensional (airfoil) flutter analysis 

is to assess the validity of the three dimensional aerodynamic theory presented in the 

proceeding sections. For this reason the two dimensional airfoil is approximated 

by a thin rectangular wing with high aspect ratios . For the current problem 

rectangular wing with AR=16 can be considered as a two dimensional wing as being 

verified in Figure 30. In Figures 30a, 30b, 30c and 30d, the curves of C Lk , C u , 

Crth r anc ^ as functions of the aspect ratio for a thin rectangular wing with 

K~ 0«25 are presented. It is shown that most curves converge at AR = 16. Values of 
SY 

CtH ' Cla 1 Cm!) r and versus K varying from 0.1 to 1 .5 are obtained with X^ £ ",2C^ 
NX— 8, NY— 10, NW=20, LW =3.0C, 'C' = 0 o l% and AR = 16 are compared with the 
exact ones in Figures 34a, 34b, 34c and 34d 0 The numerical results are in good 
agreement with the exact ones for reduced frequencies higher than 0.15, while the 
results for lower reduced frequencies are not completely satisfactory . 

As mentioned before, by plotting the value of imaginary A versus reduced 
frequency K, the flutter speed can be located at a specific K where the imaginary A 
is zero. An example of this work is presented in Figure 31 which is related to a 
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rectangular wing with AR = 16, "C = . 1%, /A = 5, = - ,2c, = 0.5, 

X^ = 0.2 and r^ = 0.5. The flutter boundary occurs at K = 0.53 and A = iAJ J/uS' 

- 1.775. Combining Equations (4.35) (4.36) and (4.37) gives the nondimensional 
flutter speed, 2Up /q ~ 1 .42 which is in good agreement with the exact solution 
(Reference 24)? The flutter speed as a function of 4^/w ^ for the same wing with /A= 5, 
Xc<= 0.2, r^ = 0.5 and X^ = “•2c is presented in Figure 32a and compared with the 
solution given by two dimensional airfoil theory (Figure 9-5(A), Reference 24). Ad- 
ditional results for the same wing with 10, X ^ = 0.2, r^ = 0.5 and X^ = “*2c 
are presented in Figure 32b. 


*See Figure 9-5(A)(j) of Reference 24 (p. 540), which yields 2 Up/cu^ - 1.4 7. 
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SECTION 8 


CONCLUSIONS 


8.1 INTRODUCTION 

The general comments of the present method are made in Subsection 8.2. The 
outline of the computer program is introduced in Subsection 8.3, while the CPU time 
and the memory space required for the computer program are listed in Subsection 8.4 

8.2 GENERAL COMMENTS 

The results obtained in the previous sections are in very good agreement with 
the existing ones. Especially, compared with the exact solutions, several results are 
remarkably accurate. For example, values of Cju t evaluated for a supersonic delta 
wing in Subsection 5.4.3, is accurate up to fifth significant figure. In addition to the 
accuracy, the method is very efficient, extremely general and very easy to use. The 
results obtained include a variety of problems for either steady or oscillatory, subsonic 
or supersonic flow. Two results for wing-body configurations in subsonic and supersonic 
flow are included, while the method can be further applied to wing-body-tail combina- 
tions, or more complex configurations. A preliminary result of flutter analysis is 
presented in Section 7. The pressure is evaluated by finite difference method. 
However, a hyperboloidal distribution of the velocity potential,^ , within each element 
in terms of at corners can be obtained by evaluating <j> at each corner by taking 
average of at the centroids of the surrounding elements. Therefore, an increase of 
generality and efficiency can be attained. In addition, the present method for the 
oscillatory flow can be generalized for the study of the unsteady flow without additional 
increase in computational complexity (Reference 25). Also it should be mentioned that 
the computer code presented here is believed to be the only one available for supersonic 
oscillatory flows around non-zero-thickness configuration. 
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8.3 COMPUTER PROGRAM 

The formulations in the previous sections are imbedded into a general computer 
program called SOSSA ACTS (Steady and Oscillatory Subsonic and Supersonic Aero- 
dynamics for Complex Transportation System). The listing of this program is 
given in Appendix E. In general/ the program is separated into three parts. The first 
part of the program is to define the geometry of the aircraft, wake and diaphragm (if 
necessary) and divide their surface surfaces into a number of elements. Each element 
is then replaced by a hyperboloidal surface defined by four corner points. This part of 
the program includes Subroutines GEOMET and VEC123. The second part of the 
program is to evaluate the coefficients ChK • bhl< anc * V/jtK anc * so ' ve a s y s 1‘ em 
linear equations by using the Gaussian elimination method for unknowns-^ , or ^ . 

This part of the program includes Subroutines CEOFF and CCGELG . The third part 
of the program is to evaluate the pressure coefficient and the generalized forces from 
the computed velocity potential <j> or cj? . This part of the program includes Subroutine 

COEFPR. 

8.4 COMPUTER TIME AND MEMORY SPACE 

The following table gives the computer time used by the program SOSSA ACTS 

in terms of the total number of elements on the wing, = 4 x N x x N y . 

Computer Time (Secs) for SOSSA ACTS 

Subsonic Supersonic 


N ELEM 

Steady 

Unsteady 

Steady 

Unsteady 

4x4x4 

29 

161 

8 

21 

4x5x5 

70 

324 

15 

38 

4x6x6 

143 

543 

29 

65 

4x7x7 

268 

928 

54 

119 
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All the results are obtained by taking advantage of the symmetry (or 
antisymmetry) in y and z directions. The wake for subsonic flow has N^, = 1 for steady 
and = 30 for unsteady flow. The diaphragm for supersonic flow has = 9. The 
computer times are expressed in seconds and were obtained on an IBM 370/145 
computer of the Boston University Computer center. 

Finally, the memory space required for steady flow is approximately given by 


N words “ 9 ' 500 + N EQN + N EQN 


( 8 . 1 ) 


where N worc j s is the number of words, while is the number of equations. 

Similarly for unsteady flow, the memory space requirement is approximately given by 


N words “ 12,500+ 22 N £QN + 2N eQN 


( 8 . 2 ) 
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Figure 2a. Wing Planform 





Figure 2c, Two Types of Element— Grids for Delta Wings. 
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Figure 4. The distribution of c J4 along 2y/b> = .7 for a 
rectangular wing with AR = 1 .0, M = .2, an 
NX = NY " 7 for comparison with results of 





0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0 

(X-X )/C 

L.E. 


Figure 5. The di stribution of along 2y/b = .707 for a tapered 

swept wing with AR - 3, TR = .5, = 45°, M — .8 y =5’ 

and NX = NY = 7 for comparison with results of Ref. 13 











1.0 


T 


T 


T 


-A- NASA TN-D-344 
— Q — Present Method 






0 .2 .4 .6 .8 1.0 

x/c 

Figure 8a. The distribution of the lift coefficient, c^ , on a 
symmetric rectangular wing with AR = 3 ,r= 0.05 
tt= 5°, M = .24 and NX=NY=7 for comparison 
with results of Ref. 15. 
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ie distribution of Cp on the upper and lower 
rfaces of a symmetric rectangular wing with 
R= 3, t= 0.05,^ = 5°, M = .24 and NX=NY=6 
r comparison with results of Ref. 15. 
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Figure 10. Lift Distribution Coefficient, Ci,, Versus Aspect Ratio AR, for Delta Wing 

u,-?L e ,° d r y - Sub J on ' ,c FI ° W ' W5fh t = 0.001 , M = 0, N = N = 7. Comparison 
With Lifting Surface Theory of Reference 16. ^ 
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Figure 11a. The distribution of section lift coefficient, C„ , 
for a wing-body configuration with = 6°, 
ei B = 0°, t= 9%, M=0, b=6c, r = 0.5c and 
200 elements on the whole wing for the com- 
parison with results of Ref. 17 



Figure lib. The distribution of along three circumferential 

stations for a wing-body configuration with ^ w = 6° 
o(r>= 0°, C= 9%, M=0, b=6c, r=0.5c 
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Figure 13b, Analysis of Convergence: Lift Distribution Coefficient, c^, Versus x/c, 
at 2y/b = 0. 1328, for Rectangular Wing With Biconvex Section, 
Oscillating in Bending Mode in Subsonic Flow> for AR = 3, t = 0.01 , 

M = 0.24, & - 0°, K = 0.47, N w - 30, L w = 3.5c, NX = NY = 5,6,7. 




Figure 14. Analysis of Convergence: 
Versus x/c, at 2y/b = 0.1 
Biconvex Section, Osciilc 
Flow, for AR =3, t = 0.0 
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Figure 15. A. s of Convergence: Lift Dislribufion Coefficient, c^. Versus 
x/< 2y/b — 0.1328, for Rectangular Wing With Biconvex Section 
O:-- ing in Bending Mode in Subsonic Flow, for AR = 3, T = 0.01 

M = 4, a = 0°, K = 0.47, Ax w = 0. 1, N = N = 7. 
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Li ft - Coefficient, Ci, Versus M, for Rectangular Wing Oscillating in 
Pitch, With AR-2, t = 0.001, k= 1, N x = 10, N y = 6, N w = 20, 
Lyy/c — 2, Nq = 30. Comparison With Results of Reference IS. 







Figure/ 76. 


Moment Coefficient, C/yi, Versus M, for Rectangular Wing Oscillating 
in Pitch, for AR -2, T = 0.001, k =1 , N x = TO, Ny = 6, Ny ^ = 20, 
L\a// c = 2, Nr. = 30. Comparison With Results of Reference 18. 






18a. Thickness effects in oscillatory subsonic and 
supersonic flows. Results are total lift 
coefficient/ Cl, versus M, for a rectangular 
wing with AR = 2 and K - tuc/2U«^ — » 
oscillating in pitch about axis X - X,^. 
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Figure 18b. Thickness Effect in Oscillatory Subsonjc and Supersonic Flows. 

Results are total moment coefficient, C M , versus M, for a 
rectangular wing with AR =2 and k = 00 c/^Uqq = 1 .0 oscillating 
in pitch about axis x = x . 




Figure 13a. Lift Coefficient, Qu Versus k, for Rectangular Wing Oscillating m 
9 Piur.gc, With AR = 2, t = 0.001 , M = 0, N x = lO, N y - 6, N w - . 

•* ' « Comparison With Results of Reference 


Uk, - 2c. 
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Total lift coefficient, C|_, versus M, for a rectangular wing 
with AR = 2, r — 0. and k — ^c/21}^ = 1.0 oscil lating in 
plunge. Results are related to both subsonic and supersonic 
flows. 





Figure 20b. Total moment coefficient, C M , versus M, for a rectangular 
wing with AR — 2, t = 0.1% and k = wc/ZUqq oscillating in 
plunge. Results are related to both subsonic and supersonic 
flows (same problem as Figure 20a). 








Figure 22a. The l ift distribution on symmetric rectangular wing 

with AR = 3, T - 5%, a = 5°, M = 1 .3 and NX = NY = 7 
for the comparison with results of Ref. 15. 












Figure 23. Analysis of Convergence: Potential Distribution, $, Versus x/c, at y - 0, 
for Rectangular Wing With Biconvex Section, in Steady Supersonic Flow, 
for AR = 3, T = 0.05, M = 1.3, e = 0°, Np = 3N x , NX = NY =5,6,7. 







Figure 25. Lift Distribution Coefficient, Cja-f for Delta Wing With Subsonic 
Leading Edge, in Steady Supersonic Flow, With B/tan A = 0.833, 
T = 0, N = N =7. Comparison With Exact Conical-Flow 
Solution, X Refei^ence 20 . 



Figure 26. Study of existence of diaphragm; the distribution of velocity potential, 
0, for a circular cone with subsonic leading edge, r = 1 , 1-5, 

M = 2.0 and NX = 12, NY - 5. Results obtained with and without 
diaphragm are compared. 
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Figure 28. Lift coefficient,^, for a rectangu 
k = tuc/2U 0O = 0.1, M = 1 .3, AR 1 
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29 Analysis of Convergence: Distribution of? =5e l Versus x/c, at 
2y/b — 0*5, for Rectangular Wing With Biconvex Section, Oscillating 
in Bending Mode in Supersonic Flow, for AR = 3, T “ 0.01, M - 1 .3, 
K = 0. 1 , a = 0°, N n = 3N , NX = NY = 5,6,7. 
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' Figu# 30c. as a function of aspect ratio for a rectangular wing with 

. T = 0.1% / M = 0 and NX = 8, NY = 10 (X £A = -0.2C). 
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Figure 30d . 
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Figure 31. Parameter A = oo^/a)^ as a function of reduced frequency, K - ® c/2\J^ f 
for a rectangular airfoil with AR = 1 6, M = 0, = 0.5, X a — 0.2, 

= 0.5, M = 5 and NX = 8, NY = 10 (X = -0.2C). 



T 


two dimensional airfoil 


© present method 


Figure 32a. Flutter speed as a function of ^ for a rectangular wing with 
AR = 16, M = 0, t = 0. 1%, p, = 5, Xa = 0.2, = 0.5, and 

NX = 8, NY = 10. Results are compared with exact solution 
given by two dimensional airfoil theory (Ref. 24) (X c * = -0.2C 


two dimensional 


q present method 


Figure 32b. Flutter speed as a function of (Ul/oJq, for a rectangular wing with 
AR = 16, M = 0, T = 0.1%, u= 10, Xa = 0.2, = 0.5 and 

NX = 8, NY = 10. Results are compered with exact solution 
given by two dimensional airfoil theory (Ref. 24)(Xr* = -0.2C). 
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Fiaure 33a. Cl h as a function of the reduced frequency, k ~® c / 
rectangular wing with AR = 16, t = 0.1%, M - 0 an< 
NY = 10. Comparison with the exact solution given 
theory (Ref. 24) (X £A = -0.2C). 
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Figure 33b. Cl^ as a function of the reduced frequency, k = U)c/2Uoo for c 
rectangular wing with AR = 16 , r = 0. 1%, M - 0 and NX = 8 
NY = 10. Comparison with the exact solution given by 2-D 
airfoil theory (ref. 24) (X CA = -0.2C). 
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Figure 33c. C M |i as a function of reduced frequency, k =<W c/^Uqo for a 

rectangular wing with AR = 16, t = 0.1%, M = 0 and NX = 8, 
NY = 10. Comparison with the exact solution given by 2-D 
airfoil theory (Ref. 24) (X^ = -0.2C). 






2.0 


air run incur/ 



Figure 33d. C tAa as a function of reduced frequency, k = Wc/IUqq for a rectangular 

wing with AR = 16, t = 0.1%, M = 0 and NX = 8, NY = 10. Comparison 
with the exact solution given by 2-D airfoil theory (Ref. 24) (X p . = -0.2C). 
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APPENDIX A 

FORMULATIONS FOR STEADY SUBSONIC FLOW 
A.l INTRODUCTION 

The derivation of Equations (3.10) and (3.11) is presented in Reference 6. 

In this appendix it will be shown by differentiation that doublet and source integrals in 
Equations (3.10) and (3.11) are valid for any planar quadrilateral element. In Sub- 
section A. 2 f basic equations for hyperboloidal surface are introduced. In Subsection 
A. 3, the doublet integral in Equation (3.10) is shown to be valid for any planar quadri- 
lateral element. .In Subsection A. 4, the source integral in Equation (3.11) is shown 
to be valid for any planar quadrilateral element. In addition, tne wake prm in 
Equation (3.16) is derived in Subsection A. 5. 

A. 2 BASIC EQUATION 


According to Equation (3.15) 
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For hyperboloidal element’s with aj and a, 2 defined by Equations (3.12), (3.13), 
and (3.T4) with 

A % - ) a cu | ^ J 

the doublet and source integrals in Equations (3.4) and (3.5) can be rewritten as 
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since 
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A. 3 DOUBLET INTEGRAL 

The following proves that Equation (3.10) is valid for any quadrilateral 
element. By using Equations (A.3) to (A. 7) one obtains: 


fth' 


-/ ~C % 1 — 


I 


tan __ _ 

( ( + f (f **.)•» *«J f 

W‘f Ci *Zx «.») 


A- 2 






r~H** at)* 


— ., _ _ _ _ ^ x ., _'_ ■ — ^' ( ^J 

t| • t 5 X * 2 ^} 

= 7 . .... '•,»., > - y [((^•(|x«-)(j.|Kp x 

( ^ • X, * <T* ) + (c <£x p 3 ) • (| * «.2)( • <y * *) ” 1 *• 1 * * ^ ^ V i ^ ^ t * ' 6 x ft $ $ ’ ? | ^ A ’ 1 2 ^ 

Considering terms inside the bracket yields 

&>* a i) * <u)^ # ft) ” ^ ^Hft 11 ^ ; . 

0-iK^* ft *Kft x a,)) ( J * J ) 

*|(CvtKV**) ' C a*- aOU.* f )} 4 • V-Ccf- j )t v lO-^Kft* 

liu-f )] l V' «> 4 *> + (< V \) 1 ft • ■ ( j- <*) *• v f)> ) 

ifr *■**»)' Cn-V ( «.-^- c V^> ( V 



(c <U‘ fj) <-«$* $• *«,xij( j- 

since (note the change of the order in the triple product) 

-(V^vf 3 )cr^^-c- v 
= oT. f )&<»'*>)<$ • ^x^,)'Ca 2 '^) C a 2 * |x a,)) 

-Ccf-^)(«.'a.)-c| •ctO*J C^' fVv 


_ _ , _ — 


1 y H Y ft* ^ 


(A. 14) 


Finally, by using Equation (D.2) and combining Equations (A. 12) and (A. 13) one 
obtains 


-l 


4 S 





J=L -k 

4» 11*1.1 If* i| 



CiiOiiiM ^ |X ~v i 

j v^r 


(A. 15) 


A - 4 



Finally, by using Equations (A. 3), (A. 7) and 
'd . . - — 


- o 


Ti^Y at) ' a,/a ' 

one obtains. 


’i”> 7 i*l P o/fTf-v-J, 


_ 2 L r 

*5 ^vt- ly**i* 


*v t 

( c V^ >C T ( f * p c ^ ^ ^ x ^ 




(A. 16) 


(A. 18) 


s ft * fli* «■» 

. ( vi )% 

This proves that Equation (3.10) provides an analytical solutions of the doublet integral, 
given by Equation (A. 9), for any hyperboloidal element. 


A. 4 SOURCE INTEGRAL 


Use Equation (A. 3) and note that 
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and combining Equation (A. 14) with Equations (A. 5) and (A. 6) yields 
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Finally combining Equations (A. 21), (A. 22) and (A. 24) yields 
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This proves that Equation (3.11) provides an analytical solution of source integral. 


given by Equation (A. 10), for any quadrilateral planar element. 


A. 5 WAKE COEFFICIENT 

Consider Figure 2b, and assume that the wake element is truncated such that 
°1 = Pi = = (A. 26) 

By letting^ go to infinity, wake coefficient can be obtained from Ip in Equation 
(3.10). Note that (see Figure 2b) 
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If is convenient fo separate the contribution from the trailing edge ( J 
and the edge that goes to infinity =1): 
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APPENDIX B 


FORMULATIONS FOR STEADY SUPERSONIC FLOW 
B.l INTRODUCTION 

The derivation of Equations (5.8) and (5.9) is presented in Reference 7. 
In Subsections B.2 and B.3 it will be shown that these equations are valid for any 
quadrilateral planar element. In Subsection B.4, formulations of the doublet and 
source integral of the element intersected with the Mach forecone are considered. 


For hyperboloidal elements with q, aj, cu, and cte defined by Equations 
(3.12), (3.13), and (3.14) and (A. 8), the doublet and source integral in Equations 
(4.7a) and (4.7b) can be rewritten as 
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DOUBLET INTEGRAL 


Following is the proof that Equation (5.8) is valid for any quadrilateral 
element within the Mach forecone. By using Equations (A. 3), (A. 4) and 
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Next, consider the second mixed derivative, since 
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Thus, prove that Equation (5.8) provides an analytical solution of the doublet integral 
given by Equation (B.1) for any quadrilateral element. 


B.3 SOURCE INTEGRAL 


Following is the proof that Equation (5. 19) is valid for any quadrilateral 
element within the Mach forecone. Noting that,, as shown in Appendix D, 
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In addition, from Equation (B.9) 
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Finally, combining Equations (4.10), (B.12b), (B.13) and (B. 14) and noting that 
10 <1 = -nor i yields 


? a I-6 J_ > 

Tf J)0 


TT h^h* ( ^ 11^(1 || y| 






K" s 7T TT [ 3 ia.*«» l 

According to Equation (D,7) for 


J X^i-H & X 


■ I aTj( OU. — ^<7X ^i. 61 *>) n(^-q ( xaj^) 






flO + { V a » x< *^ 

= * ,x «*^ +t V a,x 

- ^ n 0 n I Q| X o.J\ 

By using Equation (B.16), Equation (B.15) reduces to 
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This proves that Equation (5.9) provides an analytical solution of the integral 
in Equation (5.5) for any quadrilateral element. 

B.4 FINITE PARTS OF INTEGRALS 

In order to extend the results to a planar quadrilateral element intersected with 
the Mach forecone, the finite part of the following special integral is investigated. 
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where g is a regular function and H is the Heaviside function. Equation (B. 18) can be 
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This shows that the singular contribution disappears and should not be taken into 
account. * 

Thus, consider the source integral given in Equation (B.9). According to 
Equations (B.15a) and (B.17), the integran of Equation (B.9), , can be 
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First, according to Equation (B . 12b), Equation (B.22b) can be rewritten as 
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Note that'v/^o^ can be expressed as ' where l\(/j ) is a regular 

function and is defined such that ^ 0 ’^ • Compared with Equation (B.18), it 
can be concluded that the portions along the intersection line of the element with the 
Mach forecone yield no contributions to the integral S, • Therefore, in Equation 
(B.22a) 
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if the .edge of /j*t or ^=-| , respectively, is completely outside the Mach forecone. 
Otherwise, I S1 is given by Equation (5.10) if the corner point is inside the Mach 
forecone, or 
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if the corner point is outside the Mach forecone. Similarly, Equation (B.23b) can be 
rewritten as 
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if the edge of-| - j , or ~^s-\ , respectively/ is completely outside the Mach forecone. 
Otherwise, is given by Equation (5.11) if the corner point is inside the Mach 

forecone, or 
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if the corner point is outside the Mach forecone. Finally, considering Equation 
(D. 10), Equations (B.24b) and (B.7) can be rewritten as 
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if the edge of 1 or , respectively, is completely outside the Mach fore- 

cone. Otherwise, 1^ and Ip are given by Equations (5.12) and (5.8), if the corner 
point is inside the Mach forecone, or 

*7*) - ( yyi ) (B.39) 
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Note that the values of ^*and *7* in Equations (B.32), (B.28), (B.39) and (B.40) are 
evaluated such that =o . 

B.5 INTEGRAL EQUATIONS FOR DIAPHRAGM-ATTACHED CONFIGURATION 

As mentioned in Subsection 5.3, if diaphragms are used, both values of the 
velocity potential and its normal derivative, and ^ , of the diaphragm element 
are unknowns while two integral equations are obtained for each diaphragm element. 
Therefore, a system of equations, instead of Equation (5.3), may be obtained as 
follows. 
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centroid of the element on the upper (or lower) surface of the aircraft, while lf? p 
and are those quantities for the diaphragm element , and J ^ is given by 
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and kf^are given by 



and 
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where the superscript u (or i ) reminds that both the control point and dummy 
element are on the upper (or lower) surface and first subscript a (ord) reminds 
that the control point is on the surface of aircraft (or diaphragm) second sub- 
script m (or n) reminds that the dummy element is on the surface of aircraft (or 
diaphragm). For example, C U represents the submatrix formed by the doublet 
integrals, given by Equation (5.4), with control point on the upper surface of aircraft. 
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and dummy element on the upper surface of diaphragm, represents the submatrix 
formed by the source integrals,, given by Equation (5.5) with control point on the lower 
surface of diaphragm and dummy element on the lower surface of aircraft. 

Furthermore, if the problem is symmetric, the normal derivative of the 
velocity potential , ^ , on the diaphragm is zero, while the velocity potential, ^ , 
is unknown, 
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For the case that the problem is antisymmetric, , is zero while is unknown. 
Then Equation (B.41) reduced to Equation (B.46) with [b^l given by (B .48) and 


[°hk^ given by 
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APPENDIX C 


C.l 


FORMULATIONS FOR OSCILLATORY FLOW 
BOUNDARY CONDITION 

Tfie boundary condition for oscillatory flow is given by 

^ 5 
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0.1.1 Subsonic Flow 


By introducing variables X, Y, Z f T/-42-and<£> given by Equations (2.2) 
and (3.3)/ and noting that^-l-M^ Equation (C.l) yields 
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Small terms can be dropped out by investigating the order of each quantity. Assuming* 
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Also, [see, for instance, Equations (C.16)and (C.17)], it is noted that 
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Investigating the order of each term of Equation (C.5) by Equations (C.6) to (C.13) 
and ignoring the terms which contain e'^^ (of order £.^), one obtains the stationary 


and the oscillatory boundary conditions which are given as 




M2 2h. 


'Z'X 


0 


(C.14) 


and 

VS- • + v n*“s- 'Vh + Tf 

^ ^ 'yX <dt vx ‘It ' 


(C. 15) 


2 . 3 

Neglecting terms of order ^ in Equation (C.14) and terms of order £ in Equation 

(C.15), one obtains 
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By using Equation (2.31), Equation (C. 22) yields 
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Equation (C.20) gives the value of to be used in Equation (4.3) 
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C . 1 . 2 Supersonic FI ow 

By introducing variable x, y, z, T,-^-and <\> given by Equations (2.2) and 
(3.3), and noting that £ 2 = M 2 - 1/ Equation (C.l) yields 

(C.25) 


_ c . v* . \ 21.+ S4. 

V xrsS’ V xr2 T M a-J* P -aX p* ^x 


with 57* defined by 
xyz 


%n 


•?x 


* + —v I +• — ~ "K 

4 J -a 2 


(C.26) 


o' 


o' 

l 

/ 


By replacing V xyz with v* yz , ^ xyz f 0 with V * yz f 0 and V xyzf with V xyz f 
Equations (C.5) to (C.22) are valid for supersonic flow. Therefore, the boundary 


condition for supersonic flow can be rewritten as 
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By using Equation (2.39), Equation (C.27) yields 
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Equation (C.28) gives the value of to be used in Equation (6.3). 
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C.2 PRESSURE COEFFICIENT 

The pressure coefficient is given by the linear Bernoulli theorem. Equation 
(2.41), as 
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C.2.1 Subsonic Flow 


By using Equation (2.31) 
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Cp for oscillatory subsonic flow can be rewritten from Equation (C.31) as 
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C.2.2 Supersonic Flow 




By using Equation (2.39) 
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APPENDIX D 
USEFUL EQUATIONS 


D.l INTRODUCTION 

In deriving the formulations presented in previous sections, several useful 
equations are applied. They are proven in this appendix. 

For the subsonic flow theory, the following two equations are used. 

2 7. 3- (C 

*-V ^ H + E I 

For the supersonic flow theory, the following five equations are used. 
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D.2 PROOF OF EQUATION (D.l) 
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4 Ax $x Py Py 4 Ay By C x Px “ ( AyCy # x P/ 4 Ax Px By Py) (D . 8) 


D.3 PROOF OF EQUATION (D.2) 


Following is the proof of Equation (D.2): 

if* «»rif * ' ( v % ) c i'^ x ^AC ( ^^)*^ K ^)} I 

= Ccf-^Xa r ap-(^vaP 2 X^'f)(Az'^-(f'^~^'^)(f'a ( x^ z ) 2 

"C ( f ‘f ) C«i ’ «») ' ^.X %' *$) 

= C V %)’Ca, •*.)<* -A.*) + ( f 

"(f t Xf- ■«iMD 2 -(f' f) 4 ( ‘ i/if ' «>) a+2 ( f- fXv^Xi'«"Kf'al; 

* 

= l V J‘ > ' £f'f v <^)( - ((«.-«.)(f *a/-i‘ 

■i (%' ^cf •a/'-2Ca J -ocf .ajtf • j[i)jJ 


D -2 


v ✓ 



(D.9) 


- o 

since, by using Equation (D.1), 

| fxU l y?z)! 2 = Cfxia l )in>)*Cf 

• «|Xa*H %-ap <K z y 

and by using ^ c) s alA-c)-C^.p) ' y Je,ds 

D.4 PROOF OF EQUATION (D.3) 

Note that 

(cu fe ) 0 ( cx 7) 

-^2 !>y)C^y 4 "Cgdy) ” CCibx~(Xx\>i)(.C2& 1 c'Cx < ^z) 

- (^x^y ~ ^Y d x ) 

> ^-Y^h Cydg + ^zbyCgciy " ^Tyl?g C 2 dy “<^^|?Y^Y^ 2 

- ^2 by £? d x “ b? d 2 + ^2 by dx + $.x bg Czdx 
' a x by c x d y - fty t> x CycJ x + 0. x tyCy d x + Ay b x C x dy 

while _____ _____ 

(a.. 0 t )C b o 3 ) -Ca. e> d K b © C > 

r ( a x c x ' ^y c y -^UCgXbxd* - by dy - bg ) 

^CAydx ' ft ydy -^ 2 d£ )(bxCx -by Cy - b 2 ^2 ) 

— O-iCi C bx dx — by d y - ) —ft* dx ( b^x “by Cy — b 2 Cz } -f* 

D - 3 


= a*c x (y 4 ^Mr'Mz>-%Cr(^ J *"^X~ l ’ i! ‘ j2 - ) ' 


(D.10) 
(D .11) 

(D. 12) 



+ flydy (bxCx-k c Y -Va C*) tCUdg Cb* Cx-byCy- t> 2 C 2 > 

- (LybiCy^z + ftx kCjdY - %UCzdy-^2oy c Y^2 —( ??:ix c 2dx _ ^^^^2 

(Q , ] 

+ fl.2bx^d2 ~ ft* by Cyciy - &y bx Cyd/ + Ct^by Of dx +ft.ybxCxdy +^b2^Hdx 
D .5 PROOF OF EQUATION (D.4) 

_ _____ t. __ __ n 3- 

r ( f® t ^ C f° t ) ( V 0 V v f )- ( ®<Tx*pJ 


+(<*•* >U,**0 0 Av)] 2 ’ 

' *'° «*■) u ) 2 ) “ ( %* \ )(f 3 ?*>((& «7>. [Qi o ~%y^° 

%) -(fl2 o] +(ft®^A<b ©J*/” 

' 2 C<^ 0^)Cq, o a2.)( oa>-)C <g.© ®a0 

: 6| r |)^VA')'if^7J Cc^®%)(^© «*;-^® «*/>! 

" <5^x a,) =^x «, © ^xa,> (|X^ ® $MO 


D.6 PROOF OF EQUATIONS (D.5) and (D.6) 

Proof of Equation (D.6) is shown below. Proof of Equation (D.5) is similar; 
therefore, it is omitted. There are three different cases of Fj. Considering the first 
case aj © 02 ? 0/ and noting that 


J_ 

*7 


II cull 


n O 


(D.14) 


and 
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= 2 (S.x’te) 0 (f x = ° 


?F» 2 f'' I „ II til I'M + <*> > 

*7 ' *1 1 «*.< h Ilf, in 1 J 


* ifeTiC^ h l»tH'^ l + r Zj*»\y*^ 


\ / ® j — _ N 

i + i®^ * tTzoa. r> 0 fe +l? 5 0<? J 


JUjII ll^ll HcuU + ^Uj. \ 


i 


)|^7.ll ill'll ll^>H 
I 

"t" 


(?«>#» ik,n +ll fll 11 «>0 '7s, 


llfll 


Consider- the second case, aj © aj = 0 

*F 9 * ( — li-tiJ 'T 

? *1 *7 L f 0 4> 

I 


— — “ II %\\ 

f 0 *> *0 6 


I 


e k 9 


% 0 <i% 1 

11 1 11 


since 


? - 


^ ( | t )®av" dv' 5 ’^' 

Consider the third case, a^ ® °2 ^ 0 


(D.15) 


(0.16) 


(0.17) 


(0.18) 
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2h z-Z-C-— St>-'(-*^-)h-4r 

?r\ ?*) L ii«JI |||mJ| } 


__ J I 4a ® 4a 

" w VT- (4^-f ii i »< «“ ii 

4*ll' P 

_ 1 IK x 7>U jjjhlL 

l^©7z) a II 

11 aV 11 _ 1 

"J ^ f *1 

Combining (B.5), (B.6), and (B.7), yields 
- 

n 

Similarly/ 

•3 p) _ 

D.7 PROOF OF EQUATION D.7 

Consider the regular vector algebra rule 
& * (I* C) =T(a;*C) ' c(7-7J 

and note 


'V 


/> On ^ O 




i 0 Q i 


(D.19) 


(D.20) 


(0 .21) 


(D.22) 


a. e ~E * t - a • J c 



One obtains 


fl c x(F’«c)= b U c ' c) - cUVEJ 
* I (a> C) - c (a®T) 

On the other hand, according to Equations (D.3) and (D.22) 

0i c * (F * t) ® Of x CB <C) 

_ of © cL fc ( TT* "c) ° ("E <c) _ (X 0 ^ b * c)} 

r cl o cl (b o (X x C) - (a. • b * cJ 

while, according to Equations (D.3) and (D.22) 

( b ( SC o Z ) - rc ZC o*B) 3 ® C&l «> T) " r C 57 G 

= ^ 33 G c a>b + ZoCCHolf 

- aoc cixc’Hx^as'B (cxlerx^) 

Combining Equations (D.22), (D.23) and (D.24) yields 
a ° a. (¥xt?° “ ( a *“5 x'O 1 ' 

= (a c * tX °( a c * xTi)l 
= (% (c[ 0 t) - rca ®*E;J ® (S ~ c (G. * b>) 

- a o C o~Ex^C> •+• rcx“B o’C’xa.^ 


(D.23) 


(D.24) 


(D.25) 
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APPENDIX E 

LIST OF COMPUTER PROGRAM SOSSA ACTS 
MAIN PROGRAM 


COMPLEX AA, B3, SOURCA, SOURCB » VELPOT 

COMMON/Z ZZ1/MX ( 20 ) , NY ( 20 ) , NX Y ( 20 ) , NW ♦ KSYMMY , KSYMMZ YNSYMMY , NSYM” * 
C0MM0N/ZZZ2/T AU , SPAN, TANGLE » TANGTE , CHORD » I ZZZ* UMACH,REFL EN 
C0MM0N/ZZZ99/KPRT NT (10) ,NREAD,NWRITE 
C 

DIMENSION YK ( 3, 16, 16, 5 ) » PC ( 3 » 100) » YPP ( 3 » 100 ) ,YPM{3,100 )tYMP(., ’ 1 
l) ,YMM(3» 100) ,KWAKEl 100) , SOURCA < 100 ), SOURCB ( 100 ), AA { 10000 ) ,V = LP-: 7 
1100 ) » BB( 10000 ) 

C 

EQUIVALENCE ( YK ( 1, 1 , 1 , 1 ) , AA ( 1 ) ) 

C 

NREAD=5 
NWRITE=6 • 

NXMAX=15 

NYMAX=15 

NTMAX=100 

NXMAXP=NXMAX+1 

NYMAXP=NYMAX+1 

C 

READ! NR E AD, 2 INCASE 
2 FORMAT! IX, 12) 

DO 999 1 CASE= l , NCASE 
WRI TE (NWRITE, 1) 

1 FORMAT (• 1' /////) 

NL00P=1 

CALL COODPT ( NTOTAL, NXMAXP ,NYMAXP , YK ) 

CALL CHECK ( NTOTAL, NXMAX* NYMAX, NTMAX) 

NT2S- NTOTAL **2 «- 

CALL GEOMETt NTOTAL, NXMAXP, NYMAXP,YK, YPP, YPM, YMP, YMM.K'WAKE ) 

CALL V EC 123 ( NTOTAL » PC » YPP , YPM, YMP , YMM) 

CALL COEFFI NTOTAL, PC, YPP»YPM, YMP, YMM, K WAKE, NT2S, AA, SOURCA, SO JH C j 
150 CONTINUE i 

DO 160 NNN = 1» NT2S 
160 BB ( NNN)= AA ( NNN ) 

CALL S0LUTN(NT0TAL,NT2S,BB, SOURCA, VELPOT) 

CALL COEFPRlNTOTAL, PC, YPP, YPM, YMP, YMM, SOURCA, VELPOT) . j 

IF( IZZZ.lt. 500) GO TO 999 i 

I F.l NLOOP .EQ • 2 )G0 TO 999 { 

WRITEINWRITE, 1) I 

DO 200 1=1, NTOTAL ! 

200 SOURCA ! I ) =SOURCB( I ) i 

NL00P=2 | 

GO TO 150 1 

999 CONTINUE f 

STOP ' | 

END . ... . ... j 

' . . ■ • it 
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c 


in 


19 
10 

20 


18 

28 


12 

22 

13 

23 


15 

25 


C 


SUBROUTINE COODPT INTOTAL , NXMAXP , NYMAXP , YK) 

DOUBLE PRECISION TODAY 

COMMON/Z ZZ L/NXI 20 ) » NY ('20 ) »NXY ( 20 ) ,NW, KSYMMY » KSYMMZ « NSYMMY * NSYMMZ 
COMMON/Z ZZ2/TAU « SPAN , TANGLE , TANGTE , CHORD , I ZZ Z , UMAC H, RE FL EN 
COMMCN/ZZZIO /KWAKES ( 20 ) » NWAKE ,KXI NCR 

COMMON/ ZZZ 11/ ALFA , ALF ABC , REDFRE » R * XLEZ » XTEZ , XNOSE, XTAI L,KCM 
C 0 MMGN / Z ZZ 1 2 / NS FX , N S BO D Y , N S , NT 1 20 ) , KMORML ( 20) , KQ1AFC20 j , I S FACE ( 20 ) 
C0MMQN/ZZZ13/DCHQRDI 10 ) * CHAXI S ( 10) 

C 0 M MO N / Z Z Z9 9 / K P RI NT l 1 0 ) , N R E A D , N WR I T E 

DI MENS ION ACHORD (20), WNP P ( 2 ) , WNP M ( 2 ) , WNMP ( 2) , WNMM l 2 ) 

DIMENSION Y K ( 3 , NX M AXP » NY MAX P , 5 ) , GG ( 20 ) , KS F AC E ( 20 ) 

DIMENSION DFPP (2) ,DF°M( 2 ) ,DFMP (2) , DFMM ( 2 ) 


NSFX=12 

CALL DAT E( TODAY ) 

WRITEINWRITE, 111) TODAY 

FORMAT ( / //IX , ' DATE : ',A8//) 

DO 19 1=1,3 
READl NREAD, 10 ) GG 
WRITEINWRITE, 20) GG 
FORMAT (20A4) 

FORMAT (IX, 20A4 ) 

READINREAD, 18 ) GG( 1) ,GG( 2) , NS , KSYMMY, K SYMMZ 
WRITEtN WRITE, 28)GG( 1) ,GG(2) , NS , KSYMMY , KSYMMZ 

READ( NREAD, 18 ) GGl 1) , GG ( 2) , NW ING , NB0DY1 , NB0DY2 , N30DY2 , NDI AF 1 
1,NDIAF2» NDIAF3 

WRITE { NWR ITE , 28 )GG( 1 ) , GG ( 2 ) , NWING ♦ NBODY 1 , NB0DY2 , NB0DY3 , NDI AF1 
1,NDIAF2, NDIAF3 
FORMAT ( 2A4 , 1015) 

FORMAT { 1X.2A4, 1015) 

READINREAD, 12) GGI 1 ) , GG I 2 ) , UMACH, REDFRE 
WRITEINWRITE, 22)GG( 1) ,GGI2) , UMACH, REDFRE 
READ! NREAD, 13 ) GGI 1) , GGI 2) , ALFA, ALFA3C , I ZZZ , KCM , MZLOOP 
WRITEINWRITE, 23) GGI 1 ) , GG ( 2 ) , ALFA, ALF ABC , I ZZZ , KC M, MZLOOP 
FORMAT (2A4,3F8»3,3E12«5) 

FORMAT ( IX , 2A4 ,3 F8.3 , 3E12 • 5) 

FORMAT 1 2A4, 2F8 »3, 31 5) 

FORMAT! 1X,2A4,2F8.3 ,315) 

READINREAD, 12) GGI 1) ,GG( 2) , SPAN, XLEZ, XTEZ, TANGLE, TANGTE ,TAU 
WRITEINWRITE, 22 )GG ID ,GG( 2) , SPAN, XL EZ , XTEZ ,T ANGLE , TANGTE , TAU 
READINREAD, 12) GGI 1 ), GG ( 2 ), XNOSE , XTAI L , R 
WRITEINWRITE, 22) GGI 1) , GG I 2) , XNOSE , XTAI L , R 
FORMAT I 2A4, 8 F8 « 3 ) 

FORMAT 1 1 X , 2A4 , 8 F8 • 3 ) 

IFIUMACH.GT.I. .OR. REDFRE. EQ.O. OR. NDI AFi.NE.O*OR,NDlAF3,NE.O) 
1GQ TO 30 

READINREAD, 18 ) GGI 1) ,GG(2) , NWAKE , KXI NCR 
WRITEINWRITE, 28) GGI 1) ,GG(2) , NWAKE , KX I NCR 
CONTINUE 

READINREAD, 18) GGI 1 ) ,GGI 2) , ( KPRI NT I K ) ,K=1 , 10) 

WRITE IN WRITE, 28)GG( 1) » GG I 2) , I KPRI NT I K) , K=1 , 10 ) 


DO 21 I S F IX= 1 , NSFX 
K5FACEI ISFIX )=0 

• E“2 



21 


o o o o 


REFLEN=l • 

Ci!ORD=XT EZ-XL E Z ■ 
8ETA=SQRT(ABS(U MACH** 2-1.) ) 
KS = 0 

NSBQDY=0 

NSBTOT=0 


WING 


C 

14 

C 

11.7 


1131 

1132 

1133 


I F(NW ING.EQ.OJGO TO 199 
DO 198 T WING=1 tNWING 
KS=KS + 1 

NSB0DY=NS80DY+1 

ISFIX=0+IWING 

READ ( NREAD » 18 ) GGC 1 ) , GG ( 2 ) , NX ( KS ) »NY( KS ) » KNQRML ( KS )» KD I AF ( KS ) 
1 KWNELE t KWNSHP » KWNTYP 

WRITE ( NWRITE , 28 ) GG( 1) » GG( 2 ) ,NX( KS ) ,NY(KS) , KNQRML ( K S) ,KDI AF(KS) 
1 KWNELE, KWNSHP, KWNTYP 

NXY ( KS ) = NX ( KS ) *NY ( K S) 

KSFACEI ISF1X )=KS 
ISFACEtKS ) = I SF I X 
NT( KS )=NSBTOT 
NSBTOT=NSBTGT+NXY(KS) 

IF ( KWNTYP. EQ. 1 )GO TO 14 

READINREAD, 15) GG( 1) , GG ( 2 ) » WNPP , WNPiM , WNMP , WNMM 
WRI ?E (NWRI TE , 25) GG ( 1 ) , GG ( 2 ) , WNPP , WNPM , WNMP , WNMM 
CONTINUE 

I F ( I S FIX • EQ • 1 ) S IGNZ-+ 1. 

I F ( ISFIX • EQ. 2 ) S I GNZ=- 1 . 

WRITE (NWRITE, 117 ) ISFIX, XLEZ ,XTEZ 

FORMAT l //2X, ' PART « , I 2 , 5X , • X I N I T= « , F6. 1 ,5X , • XF I N = « , F6. 1 ) 

Ml NGLN=S P AN— 2 •* R 
WlNGWD=WINGLN/2. 

RLE=R 
RTE = R 

DXX=1 • /NX ( KS ) 

DYY=1./NY(K$) 

NXP=NX ( KS ) +1 
NYP=NY ( KS ) +1 
DO 118 I X= 1 , NXP 
DO 118 IY=1,NYP 

XX=( IX- 1 )*DXX , 

YY - ( I Y— 1 j*DYY 

GO TO (1131,1132, 1133, 1134) , KWNELE 
CONTINUE 
CSJ=XX 
eta=yy 
GO TO 113 
CONTINUE 
CSI=XX*XX 
ETA=1.— (l.-YY) **2 
GO TO 113 
CONTINUE 

f S T = X X 

ET A= 1 .- ( 1 . - Y Y ) * * 2 
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GO TO 113 
1134 CONTINUE 
CSI=XX*XX 
ETA=YY 

113 CONTINUE 

C I F { KWNTYP • NE . 1 ) GO TO 1139 

Y=W INGWD^ETA+R 
XLE=XLEZ + TANGLE* (Y-RLE) 

XTE^XTEZ+TANGTE* ( Y-RTE ) 

DCHOROl IY )-XTE-XLE 

CHAXI S( IY)=XLE+DCHORD( IY)*0.5 

X- XLE+ ( XTE-XLE ) *CS I 
C 

1139 CONTINUE 

IF ( KWNTYP .NE • 2 ) GO TO 1149 • 

X I - ( WN P P ( 1 )— WNMP( 1) )*C S l + WNMP ( 1) 

Y 1= ( WNPP ( 2 ) -WNMP { 2 ) )*CS I+WNMP ( 2) 

X2 = ( WNPM ( 1 ) -WNMM ( 1) )*CSI +WN.MM ( 1 ) 

Y2=( WNPM{ 2 )-WNMM ( 2) ) *CS I +WNMM ( 2 ) 

X3= ( WNPP ( 1 ) - WNP Vi ( I) )* ETA + WNPM ( 1) 

Y3=(WNPP( 2)tWNPM(2) )*ETA+WNPM { 2 ) 

X4=(WNMP( 1 >-WNMM{ 1) )*ETA+WNMM( 1) 

Y4= ( WNMP ( 2 ) — WNMMI 2 J )*ETA+WNMM( 2) 

OET=( Xl-X2)*( Y3-Y4)-< X3-X4)*( Y1-Y2J 
IF(DET .EQ.O. )GO TO 1147 

X«( (Xl-X2)*(X4*Y3-X3*Y4)-( X3-X4)*( X2*Y1-X1*Y2) ) /DET 
Y= ( ( Y 1— Y 2 ) * ( X 4* Y3-X3* Y4 ) — ( Y3-Y4 ) * ( X2#Y 1-X1*Y 2 ) ) /DET 
GO TO 1149 
1147 CONTINUE 
X=X4 
Y=Y4 

1149 CONTINUE 
C 

GO TO (1151, 1152» 1153), KWNSHP 

1151 CONTINUE 
C 

C FOR CIRCULAR BICONVEX ' *',• 

C 

XC=0.5* ( XLE+XTE ) 

PL=0.5*(XTE-XLE) 

H=TAU*PL 
AT / ’J=2.*H 

etAcr=i.-atau/span 

XXC=X-XC 

z=o. 

IFtH.EQ.O. .OR. IY.EQ.NYP.OR. IX.EQ. l.OR. IX.EQ.NXPIGO TO 115 
Z=STGNZ*(SQRT( ( ( PL*PL+H*H)/ ( 2.*H> >**2~XXC**2 )- ( PL*PL-H*H )/ { 2 .*H) ) 
IF(ETA.GE.ETACR)Z=Z*SQRT( l.-( (ETA-ETACR)*SPAN/ATAU)**2) 

GO TO 115 

1152 CONTINUE 

TAUBAR=TAU*.75*SQRT(3. )*(XTEZ-XLEZ I 

Z=S IGNZ*TAUBAR*SQRT (CSI )* ( 1 .-CSI)*SQRT ( l.-ETA**2 ) 

GO TO 115 

1153 CONTINUE 

Z=SIGNZ*4.*TAU*(XTEZ-XIEZ)*C$I*( l.-CSI )*SQRT(I.-ETA**2 ) 

115 CONTINUE 
C 

YK(1, IX, IY,KS)=X 
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YKI2, IX, IY,KS)=Y 
YK ( 3 ♦ IX, IY,KS )=Z 
118 CONTINUE 

KWAKES ( K S ) = 1 

198 CONTINUE 

199 CONTINUE 

NOSE 

I F ( NBQDY2 • EQ .0 ) GO TO 599 

00 598 I BODY 2= 1 »N80DY2 
KS=KS+1 

NS60DY=NSB0DY+1 
ISFIX=4+IB0DY2 

READ! NREAD, 18 ) GGI 1) ,GG( 2) ,NX( KS) ,NY( KS) ,KNORML(KS) ,KDI AFIKS) , 
1 KNSELE, KNSSHP,KNSTYP 

WRITEINWRITE, 28 )GG( 1) , GG ( 2) , NX ( KS ) , NY ( KS ) ,KNORML (KS) , KOI AFIKS) , 
1 KNSELE, KNSSHP.KNSTYP 

NXY(KS) = NX( K$ M'NYIKS) 

KSFACE 1 1 SF IX ) =KS 
ISFACE(KS)=ISFIX 
NT( KS ) =NSBTOT 
NSBTOT=NSBTOT+NXY(KS) 

IF( ISFIX.E0.5)SIGMZ=+1 
IF(ISFIX.EQ.o) S IGNZ=— 1 
NXP=NX ( KS ) +1 
NYP=NY ( KS ) +1 
X I NI T=XL EZ-XNQSE 
XF I N=XL EZ 

C WRITEINWRITE, 117) I SFI X, XINI T , XFIN 

OXX= 1 • / NX I K S ) 

XZERO=XLEZ 

RR=R 

SLOPE=XNOSE/R 
C0NST=1 • 5 

TINCR=0.5*B. 14159/NY (KS) 

DO 550 J X= 1 , NXP 

XX=(IX-1)*DXX 1 

IFIKNSELE.EQ.l >CSI=XX 

1 F I KNSELE. EQ • 2.)CS 1= XX#XX 
X = X IN I T+ 1 XFI N-X IN I T )*CSI 
IX=JX 

1 F ( KNSELE .NE • 3 ) GO TO 540 
I X =NXP- JX+1 

DX-2,.*3ETA*CONST*RR/(BETA/SLOPE + l. I 
X=XZERO-DX 
XZER0=XZER0-DX 
I F I I X • EQ. 1 ) X--XN0 SE 
540 CONTINUE 

IFIKNSSHP.EQ. I )RR=R*SQRT(CS I) 

IFIKNSSHP.EQ. 2) RR=R* I 1.0001-1 I X+0 . 01*XN0SE ) /XNOSE ) **2 ) 

I F ( KNSSHP . EQ. 3 ) RR=R*I 1 •- ( X/XNQSE 1**2) 

IFIKNSSHP.EQ. 4) RR=R* 1 1 .-ABS I X ) /XNOSE ) 

IFIKNSSHP.EQ. 5) RR=R*( l .005-ABSIX+0 .005*XNDSE) /XNOSE) 

DO 550 I Y=1 , NYP 

THE T = I I Y- 1 )*T I NCR 

THE T A= I +0 • 5*3 .14 1 59-T H E T ) *$ I GNZ 

ETA=COSI THETA ) ORIGINAL PAGE 13 
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541 

542 

543 

544 


550 

598 

599 
C 

c 

c 


992 

C 

9909 

9910 


YY=RR*ETA 

GO TO (541,542,543) ,KNSTYP 
CONTINUE 

ZZ=RR*SI N (THETA ) 

GO TO 544 
CONTINUE 

ZZ=0.01*RR*SQRT ( l.-ETA )*4.0 
GO TO 544 
CONTINUE 

zz=o* 

CONTINUE 

YK ( 1, I X » I Y » KS ) = X 
YK(2, IX, IY,KS)=YY 
YK(3, IX, I Y , KS ) =ZZ 
CONTINUE 
KWAKES (KS )=0 
CONTINUE 
CONTINUE 

— _ ---WING-DIAPHRAGM 

IF(NDIAF1.EQ.O)GO TO 999 
DO 998 IDIAF1=1,NDIAF1 
KS=KS + 1 

I S FI X=8+ ID IAF 1 
READ( NREAD, 18 ) GG(1) 

1 KDFEL 

WRI TE (NWRI TE » 28 ) GG ( 1) 

1 KDFEL 

NXY(KS)=NX(KS)*NY(KS) 

KS FAC E ( I SF I X ) = KS 
I SFACEIKS )-I SFIX 
NT (KS)=NSBTQT 
NSBT0T=NS8T0T+NXY (KS) 

DYY= 1 • /NY ( KS ) 

NYP=NY ( KS ) +1 
NXP-NX ( K S ) + 1 

XT I PC=0 • 5* ( XTEZ+WI NGWO^TANGTE+XLEZ + WI NGWD*T ANGLE ) 
DXT I P=1 • O’MXTEZ+WINGWD^TANGTE— XLEZ-WINGWD*T ANGLE ) 
XDI AC=XTI PC 
DXD I A=DXT I P 

I Ft NBDDY 2 .EQ.O ) GO TO 992 
DXD1 AF=XLEZ-XNOSE+BET A* { R+WI NGWD ) 
XTIPLE=XLEZ+WINGWD*TANGLE 
I F ( (XTI PL E-DXDI AF ) . LE .0 . ) GO TO 992 
XDI AC=0. 5* ( XTEZ +WI NGWD*TANGTE+DXD I AF ) 

DXDIA=1 .0* ( XTEZ+W INGWD*TANGTE— DXD I AF ) 

CONTINUE 

IF(NB0DY2.EQ.0.AND.NDIAF2.NE.0)G0 TO 9909 
GO TO 9910 

DXDIA=XTEZ+WINGWD*TANGTE-XLEZ-WINGWD*BETA 
XDI AC=XTEZ + WI NGWD*T AN GTE- DXD I A#0 . 5 
CONTINUE 

DIAFLN=DXD1A/(2,*BETA) 

DO 990 I Y= 1 , NYP 
YY=( IY-1)*DYY 
I F ( KOFEL E , EQ « 1 ) ET A= YY** 1 
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, GG (2) ,NX ( KS ) , NY ( KS ) , KNORML ( KS) , KDI AF(KS) , 
E , KDFSHP 

, GG ( 2 ) , NX ( KS ) ,NY( KS) , KNORML (KS) »KDI AF(KS) , 
E, KDFSHP 
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IF(KDF£LE.EQ.2)ETA=YY**2 
I F ( KDFELE • EQ .3 )ETA=YY**3 
YYY»DIAFLN*ETA 

ACHORD( IY) = ( DIAFLN-YYY)*2.*BETA*1,005 
I F ( IY.EQ.l.AND.DXDIA.EQ.DXTIP) ACHQRD ( I Y ) =DXT IP 
DXX= ACHQRD ( I Y ) /NX ( KS) • ' 

Y-YYY+WI NGWD+R 

IF( IY.EQ.NYP)Y=Y*1.002 

DO 990 I X= 1 , NXP 

X=( IX- 1 )*DXX+XDIAC-ACHORD(I Y)/2. 

I F ( IY.EQ.1.A.N0.IX*EQ.NXP)X=XTEZ+WINGWD*TANGTE 
YK(1, IX, IY,KS)=X 
YK{ 2, IX, IY,KS)=Y 
YK(3, IX, IY,KS)=0. 

990 CONTINUE 

KWAKESI K$ )=0 

998 CONTTNUE 

999 CONTINUE 


TRIANGULAR DIAPHRAGM 

IF(NDIAF3.E0.0)G0 TO 1399 
DO 1398 IDIAF3=1,NDIAF3 
KS=KS+1 

ISFIX=I2+IDIAF3 

READt NREAD, 18) GG ( 1) , GG ( 2 ) ,NX ( KS ) ♦ NY ( KS ) ,KNORML ( KS ) , KDI AF ( KS ) , 
1. KDFELE, KDFSHP 

WRITEIN WRITE, 28) GG( 1) , GG ( 2) ,NX(KS) , NY( KS ) , KNORML { KS ) ,KDI AF ( KS ) , 
1 KDFELE, KDFSHP 

NXY(KS)=NX(KS)*NY(KS) 

KSFACEI I SF I X ) =KS 
I SFAC E ( KS ) = I SF I X 
NT( KS )=NS8T0T 
NSBTOT=NSBTOT+NXY(KS) 

DYY= 1 • /MY ( KS ) 

NYP=NY(KS)+1 
NXP=NX ( K S ) 4-1 

READ (NREAD, 15) GG( 1) , GG ( 2 ) » DFPP * DFPM, DFMP , DFMM 
WRITE ( NWR ITE , 15 )GG ( 1) , GG ( 2 ), DFPP, DFPM, DFMP , DFMM 
DXX= 1 • /NX ( KS ) 

DO 1396 I X= 1 , NXP 
DO 1396 I Y= 1 , NYP 
XX=( IX-1 )*DXX 
Y Y - 1 IY-1)*DYY 
GO TO ( 1302, 1303), KDFELE 

1302 ETA=YY 
GO TO 1304 

1303 ETA=YY**2 

1304 CSI=XX. 

Xl=(DFPPC 1)-DFMP( 1) )*CS I+DFMPI 1 ) 

Y 1= ( DFPP ( 2 ) -DFMP (2) )*C$ I +OFMP ( 2 ) 

X2=( DFPM ( 1 )— OFMM( 1) )*CS I +DFMM ( 1 ) 

. Y2-( DFP'MC 2 )“DFMM{ 2) )*CSI +DFMM ( 2 ) 

X3=(DFPP( 1 )-DFPM( 1) }*ET A + DFPM ( 1 ) 

Y3= ( DFPP ( 2) -DFPM (2) )*ETA+DFPM { 2) 

X4=(DFMP< 1 ) -OFMM( 1) )*ETA+DFMM( 1) 

Y4= ( DFMP ( 2 ) — DFMM C 2 ) )*ETA + DFMM ( 2 ) 
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DET=(X1-X2)*(Y3-Y4)-(X3-X4)*(Y1-Y2) 

IF (DET.EQ.O. )GO TO 1347 

X=UX1-X2)*(X4*Y3-X3*Y4)-(X3-X4)*(X2*Y1-X1*Y2) i/DET 
Y=( (Yl-Y2)*( X4*Y3-X3*.Y4)-(Y2-Y4)*( X2*Y1-X l^Y 2 ) )/DET 
GO TO 1349 
1347 CONTINUE 
X-X4 
Y=Y4 

1349 CONTINUE 

YKd.IXf IY,KS)*X 
YK(2tIX»IY»KS)=Y 
YK ( 3 , IX, IY,KS)=0. 

1396 CONTINUE 

KWAKESI KS ) -'•/ 

1398 CONTINUE 

1399 CONTINUE 

WRITE (NWRITE, 31 ) KWNEL E , KNSELE , KBDE LE , KTLELE , KWNSHP , KNSSHP , K30SHP 
1 ,KTLSHP,KDFELE 

31 FORMAT C//5X* 'WING* , 3X, 'NOSE' ,3X, • BODY' ,3X, 'TAIL • , 3X , ' 01 APH • 74 ! 7/ 
1 417, IX, 17/) 

C 

C 

c 

IF(KS.NE.NS)CALL DE8UG(1020) 

IF (NSBOOY.LT. NS )KKD I AF=1 
I F ( KSYMMY • EQ .0 ) NSYMMY= 1 
IF(KSYMMY.NE.0)NSYMMY=2 
I F ( KS YMMZ • EQ . 0 ) NSYMMZ=1 
IF I KSYMMZ »NE .0 ) NSYMMZ = 2 

IF( KSYMMZ .NE.O .AND.KKOIAF .EQ. 1 )NSYMMZ=l 
IFtMZLOOP.EQ. DNSYMMZ-l 
C 

NTOT AL=NSBTOT 
C 

RETURN 

END 


C 

SUBROUTINE CHECK ( NTOT AL , NXM AX, NYMAX , NTMAX ) 

COMMQN/Z ZZ1/NX ( 20 ) , NY ( 20 ) , NXY ( 20 ) , MW , KSYMMY , KSYMMZ ,NSYMMY , NSYM'-'Z 
COMMON/ ZZZ 12/NS FX ,NSBODY, NS , NT( 20 ) , KNORML ( 20 ) , KDI AF120 ) , ISFACE(20) 
C 

DO 100 I S= 1 , NS 

IF (NX (IS) *GT .NXMAX) CALL OEBUG( 1001) 

I F ( NY ( I S ) • GT . NYMAX ) CALL DEBUG( 1002) 

IF (NTOTAL.GT. NTMAX) CALL DEBUG( 1003) 

100 CONTINUE 

IF(ALFA.NE.O. .AND. ( 1Z ZZ . EQ. 12 .OR . I ZZZ . EQ. 13 .OR . I ZZ Z. EQ . 14) ) 

1C ALL DEBUG (1004) 

IF ( KSYMMZ. NE.-l. AND. ( I ZZZ. EQ. 13.0R . I ZZZ . EQ. 14) ) 

1C ALL DEBUG ( 1005 ) 

RETURN 

END ORIGINAL PAGES U 

OF POOR ftUAMg, 


r: 
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c 

c 

c 


c 


c 

c 


c 

c 

c 

c 


c 

906 


SUBROUTINE 


GEOMET ( NTOT AL » NX MAXP » NYMAXP » YK» YPP » YPM » YMP » YMM,KWAKE) 


THIS SUBROUTINE IS FOR QUADRILATERAL ELEMENTS 


COMMON/ ZZZl/NXt 20 ) ♦ NY ( 20 ) ,NXY ( 20) ,NW , ^SYMMY , K SYMMZ .NSYMM Y ,NSYMMZ 
COMMON/ ZZZ2/TAIJ » SPAN* TANGLE » TANGTE * CHORD* I ZZZt UM AC H»RcFLEN 

COMMON/ ZZZ 12/NS FX,NSBODY, NS t NT ( 20) ,KN0RML(20) « KD I AF ( 20 ) » I SPACE (20 ) 
CQMM0N/ZZZ99/KPRI NT ( 10 ) * NREADtNWRI TE 4I ___ . , v 

D I MENS ION YK ( 3» NXMA XP » NYMAXP » 5 ) * YPP ( 3 ,NTOTAL) » YPMI 3»NT0T AL) « 

1YMP ( 3 » NTOT AL ) » YMM ( 3 t NTdTAL ) ,KWAKE ( NTOTAL ) 


DO 1999 I S- 1 * NS 
I SF I X= I SF ACE ( I S ) 
NXX--NXI IS) 
NYY=NY( IS) 

DO 999 I X= 1 » NXX 
DO 999 I Y= 1 ♦ NY Y 


IND=IX+NX( IS)*( IY-l )+NT US) 
IF(KNORML( IS ) .EQ.-DGO TO 906 


+- ++ 

I XMM= I X 
IXPM=IX+1 
I XPP= IX + 1 

ixmp=ix 

IYMM=IY 
I YPM= I Y 
1YPP = I Y+ 1 
I YMP= I Y+ 1 

GO TO 907 
CONTINUE 


-+ - 

++ + 


I XMM= IX 
IXMP- I X 
1XPP=IX+1 
IXPM=IX+i 
IYMM-IY+1 
I YMP= I Y 
I YPP= I Y 
I YPM- I Y+l 
907 CONTINUE 
C 

DO 908 K- 1 » 3 

YPP ( K » I NO I - YK ( K » I XPP ♦ I YPP » I S ) 

YPM ( K 1 1ND ) =YK ( K , I XPMt IYPM » I S ) 
YMP ( K » I NO ) = YK. ( K 1 1 XMP* I YMP * 1 S ) 
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YMM ( K, IND)=YK(K,IXMM, IYMM,IS) 

908 CONTINUE 
C 

KWAKEt IND)=0 

I F ( K WAKES ( IS) .EQ.l.AND. I X, EQ . NXX ) K WAKE ( IND)=1 
999 CONTINUE 
1999 CONTINUE 

I F { KSYMMY .NE • 0 ) GO TO 701 
NT0TAL=2* NTOTAL 
NS=2*NS 

NTHALF-NTOTAL/2 
NSHALF=NS/2 
DO 2990 I S=1 * NSHALF 
JS= IS+NSHALF 
KNORMLl JS)=KNORML (IS) 

KDIAFt JS )=KDIAF( IS) 

NT ( JS )=NT (IS) +NTHALF 
NXY ( J S ) =NXY (IS) 

NX (JS)=NX( IS ) 

NY ( JS )=NY ( IS ) 

ISFACEt JS)=ISFACE(IS1+100 - 

2990 CONTINUE 

DO 300 I R=l, NTHAL F 

il=ir-«-nthalf 

DO 299 K= 1 » 3 

GO TO (2991,2992, 2991) ,K 

2991 CONTINUE 
YPP(K,IL)=+YPM(K,IR) 

YMP ( K» I L )=+YMM ( K, IR) 

YPM ( K » I L ) =+Y P P ( K , 1 R ) 

YMM ( K , IL )=+YMP ( K* IR ) 

GO TO 299 

2992 CONTINUE 

YPP(K,IL)=-YPM(K, IR) ' 

YMPCK, IL)=-YMM(K, IR) 

YPM ( K » I L )=-YPP ( K» IR) 

YMM(K,IL)=-YMP(K, IR) 

299 CONTINUE 

KWAKEt IL)=KWAKEUR) 

300 CONTINUE 
701 CONTINUE 

IF ( KPRI NT ( 1 ) .EQ.L) 

1C ALL PR IN T A ( NTOTAL » NXMAXP , NYMAXP , YK » PC » YPP , YPM, YMP , YMM , 1 ) 
IF(KPRINT(2) • EQ • 1 ) 

1C ALL PRINT A { NTOTAL, NXMAXP,NYMAXP,YK f PC , YPP , YPM , YMP , YMM , 2 ) 
I F t KPRI NT ( 3 ) » EQ • 1 ) 

1C ALL PRINT A ( NTOTAL, NXMAXP, NYMAXP, YK, PC , YPP , YPM, YMP , YMM , 3 ) 
RETURN 
END 
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SUBROUTINE V EC1 23 ( NTOTAL , PC , YPP , YPM , YMP , YMM ) 

COMMON/ Z ZZ1/NXI 20 ) , MY (20) ,NXY( 20) , NW, KSYMMY, KSYMMZ »NSYMMY , NSYMMZ 
COMMON/ ZZ Z 2/T AU , SPAN, TANGLE, TANGTE , CHORD , I ZZZ , UMACH , RE FL EN 
C0MM0N/ZZZ11/ALFA,ALFABC,REDFRE,80DYR,XLEZ,XTEZ,XN0SE,XTAIL 
C 0 M M G N / Z Z Z 9 9 / K P R I N T ( 1 0 ) , N R E A D , N W R I T E 

DIMENSION YPP13, NTOTAL) . YPM { 3 , NTOTAL ), YMP ( 3 , NTOTAL ), YMM ( 3, NTOTAL ) 
DIMENSION PC (3, NTOTAL) 

C 

ALFAR=ALFA*3. 14159/180. 

SI NALF-S I N ( AL FAR ) 

COSALF=CGS(ALFAR) 

8ETA=SQRT(ABS{ l.-UMACH*UMACH) ) 

DO 100 I ND=1, NTOTAL 
YPP 1= YPP ( 1,1 NO ) 

YPP3=YPP ( 3, I ND ) 

YPP 1 1, IND)=(YPP1*C0SALF+YPP3*SINALF) /BETA 
YPP (3, IND)=-YPP1*SINALF+YPP3*C0SALF 
YPMl-YPMU, IMD) 

YPM3=YPM<3, IND) 

YPM I 1 , IND)=( YPM1*C0SALF+YPM3*SINALF ) /BETA 
YPM (3, 1ND)=-YPM1*SINALF+YPM3*C0SALF 
YMP 1= YMP ( 1 , IND ) 

YMP3=YMP(3, IND) 

YMP(1,IND)=(YMP1*C0SALF+YMP3*SINALF)/BETA 
YM P ( 3 , I N D ) =-Y M P 1* S I NA LF + YMP 3*C OS AL F 
Y MM 1= YMM ( 1 , 1 N D ) 

YMM3=YMM ( 3 , IND ) 

YMM( 1, IND)=(YMM1*C0SALF+YMM3*SINALF) /BETA 
YMM (3, IND) =— Y M M 1* S I NA LF+ YMM 3*C OSA L F 
100 CONTINUE 

DO 200 I ND=1, NTOTAL 
DO 199 K= 1 , 3 

PC ( K , IND )= l YPP ( K» IND) +Y PM { K 
PI(K,IND)=(YPP(K, IND)+YPM(K 
P2(K, IND ) = ( YP P( K» IND) -YPM IK 
P3(K,IND)=(YPP(K, IND)— YPM IK 

199 CONTINUE 

200 CONTINUE 
IF(KPRINT(4) .EQ.l) 

1C ALL PR I NTA ( NTOTAL , NXMAXP , NY MAXP , YK , PC , YPP , YPM , YMP , YMM , 4 ) 

RETURN 
END 


C 

SUBROUTINE DEBUG! K) 

C 

C0MM0N/ZZZ99/KPRINT (10) , NREAD, NWRI TE 
WR I TE ( NWR I TE , 1 ) K 

1 FORMAT (/2X, • ERROR CODE = • » 1 6/ > 

CALL EXIT 
END 




, IND)+YMP(K, INDH YMMtK, IND) 1/4. 
, IND )-YMP ( K * l ND )-YMM( K, IND) )/4. 
, INDJ+YMP (K • IND)-YMM(Kf IND) )/ 4. 
, IND)-YMP(K, IND)+YMM(K, IND) )/4. 


E-ll 


c 


c 


1 

110 


112 

113 


114 

115 


116 


2 


222 

221 

220 

3 


331 

332 
330 

4 

440 

444 

445 


SUBROUTINE PRINTA(NTGTAL»NXP,NYP,YK,PC,YPP,YPM,YMP,YMM,NPRINT) 
COM MON/ Z ZZ1/NX { 20) , NY ( 20) ,NXY( 20) , NW , KS YMMY> KS YMMZ , NSYMMY , NS YMMZ 
COMMON/Z ZZ2/T AU , SPAN , TANGLE, TANGTE, CHORD, I ZZZ , UMAC H, RE FL EN 
C0MM0N/ZZZ11/ALFA, ALFABC , REDFRE ,R , XLEZ , XTEZ , XNOSE , XTAI L 
COMMON/ ZZZ1 2/NS FX ,NSBODY , NS ,NT ( 20 ) * KNORML ( 20 ) » KD I AF { 20 ) » IS FACE ( 2C 
C0MM0N/ZZZ99/KPRINT (10) , NREAD»N WRITE 

DIMENSION YK(3, NXP,NYP,5) , PC ( 3 , NTO TAL ) , YP P ( 3 , NTOTAL ) » YPM ( 3 . NTOTAL 
1) » YMP ( 3, NTOTAL ) »YMM(3, NTOTAL) 


NY 4=4# ( NY ( 1 ) -1 ) 


GO 10(1*2,3,4 ),NPR1NT 
CONTINUE 

WRITE (NWRITE, 110) 

F0RMAT(//2X, •SPECIFICATIONS OF THE PROBLEM'/) 

DO 112 I 1 = 1, NS 
I F (NXY ( 1 1 ) .EQ .0 )G0 TO 112 
WRI TE (NWRITE, 113) I SFACEdI ) 

WRITE (NWRITE, 114) NX (ID ,NY( II ) 

CONTINUE 

FORMAT (/2X, • FOR PART ',12) 

WRITE ( N WRI TE , 1 15 ) NTOTAL , KSYMMY ,KSYMMZ, RE FLEN, SPAN , TAU , 

1ALFA, ALFAB.C, UMACH, I ZZZ 


FORMAT (2X, 'NX=» , I2/2X, ' NY= • ,12/) 

FORMAT (/2X, ' NTOTAL= • , I3//2X, ' KSYMMY® ' , 12/2X, 'KSYMMZ®' , 12// 

12X, 'REFERENCE LENGTH® ' , F6.2/2X, ' SPAN =',F6.2/ 

12X , ' WI NG THICKNESS = ' , F9 . 5//2X, ' ALFA® ' , F 7 . 3/2X , ' ALFABC® ' , F7 . 3 // 
12X , ' MACH NUMBER » • , F 7.3//2X , • I ZZZ= • , I 5// ) 

WRITE (NWRITE, 1 1 6 ) TANGLE , TANGTE , CHORD, R EDFRE , R 
FORMAT (2X, 'TANGLE® • , F6 . 2/ 2X , • TANGT E= • , F6. 2//2X , 'CHORD =• ,F6. 
12//2X, 'REDFRE®' ,F8.4/2X, » RADIUS® ' ,F6.2) 

RETURN 
CONTINUE 
DO 220 IS- 1« NS 
I SFIX® I SFACE (IS) 

NXX=MX( IS )+l 
N YY =NY ( I S ) +1 
WRITE (NWRITE, 222) 

FORMAT ( / /4X , • X • , 10X , ' Y ' , 1 OX , • Z • / / ) 

WRITE (NWRITE, 221) ( ( ( YK ( J , I X , I Y , I S ) , J = L , 3 ) , I X=1 , NXX ) , IY®1 , NYY ) 

FORMAT ( IX , F10. 5 , IX, F10. 5 , IX , F10. 5 ) 

CONTINUE 

RETURN 

CONTINUE 


WRITE (NWRITE, 332) 

DO 330 1=1, NTOTAL ' 

WRI TE (NWRITE , 3311 1 , ( YPP ( K , I ) , K=1 ,3 ) , ( YPM( K , I ) , K= 1 , 3 ) , ( YMP ( K , I ) 
1 K= 1 ,3 ) * ( YMM ( K » I ) , K=1 * 3 ) 

FORMAT ( 1 X *. 13 ,3X,3F9»5,3X, 3F9. 5,3X, 3F9. 5,3X,3E9.5) 


p ORMAT { / / 2X , 1 1 N D ' , 1 6 X » ' P P ' , 2 8 X , 

CONTINUE 

RETURN 

CONTINUE 

WRITE (NWRITE, 440) 

FORMAT (////2X, ' IND' ,4X, 'XPC' ,7X 
DO 444 1=1, NTOTAL 
WRITE (NWRITE, 445) I , (PC(K,I),K=1 
FORMAT ( IX , 1 3 , 12F10 . 5 ) 

RETURN 

END 


•PM • , 28X, 'MP * ,28X, 'MM') 


, ' YPC ' , 7 X , • Z P C ' ) 


,3) 
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SUBROUTINE COEFPR ( NTOTAL , PC , YPP ,YPM, YMP, YMM, SOURCE ,VELPOT I 
COMPLEX SOURCE, VELPOT, UN I MAG, A A, CL, CM, ACL, ACM, VEL PTE 
COMPLEX VELPTX » PHI 1 »PHI 2 * PHI 3 » COEF 1 » COEF2 *C0EF3 

COMMON/ ZZZ1 /NX ( 20) ,NY(20) ,NXY( 20) , NW , KSYMMY , KS YMMZ ,NSYMMY , NSYMMZ 
C0MM0N/ZZZ2/TAU , SPAN, TANGLE , TANGTE , CHORO , I ZZ Z , UMAC H , RE FL EN 
COMMON/ ZZZ1 1/ ALFA ,ALF ABC » REDFRE » 80QYR , XLEZ » XTEZ ,XN USE » XT AI L , KCM 
COMMON/ ZZZ12/NSFX,NSB0DY ,NS ,NT (20 ) ,KN0RML(20) ,KDIAF(20) , I S FACE ( 2C 
COMMON/ ZZZ 13/ DC HORD ( 10 ) » CHAXI S ( 10 ) 

COMMON/ ZZZ99/KPRI NT (10) , NREAD »NWRI TE 
C 

DIMENSION YPP (3, NTOTAL) , YPM( 3, NTOTAL ) , YMP (3 , NTOTAL ) , YMM( 3, NTOTAL ] 
DIMENSION SOURCE (NTOTAL) , VEL POT ( NTOTAL ) , PC ( 3 , NTOTAL ) 

DIMENSION ACL ( 10) ,ACM( 10 ) , VELPTE ( 10 ) , AA ( 1 ) 

C 

BETA2S=ABS(1.-UMACH**2) 

UNI MAG= ( 0 • , 1 • ) 

BETA=SQRT( BETA2S) 


I F ( UMACH. GT • 1 ) SGNEXP=— 1 • 

IF (UMACH. LT.1)SGNEXP=+1. 

C 

I F (KCM.EQ. 11 ) AX ISMX=XLEZ+ ( XTEZ-XLEZ ) *0 . 15 
I F ( KCM.EQ • 12 ) AXISMX=XL£Z+ (XTEZ-XLEZ ) *0 . 20 
I F( KCM.EQ. 13) AX I SMX=XLEZ+( XTEZ-XLEZ) *0.25 
IF ( KCM.EQ. 14) AX I SMX=XLEZ + ( XTEZ-XLEZ) *0.30 
. IF (KCM.EQ. 15 )AXISMX=XLEZ+ (XTEZ-XLEZ) *0.325 
I F( KCM.EQ. 16 )AXlSMX=XLEZ+( XTEZ-XLEZ )*0.350 
IF( KCM.EQ. 17) AXISMX=XLEZ+( XTEZ-XLEZ )*3. 1375 
IF (KCM.EQ. 18) AX I SMX=XLEZ+ (XTEZ-XLEZ) *0.4500 
C 

DO 600 IND=1, NTOTAL 

VELPOT( IND)=SOURCE( IN D) *CEXP ( SGNEXP*UN1MAG*REDFRE*UMACH**2 
1 *PC( 1 , I ND)/ BETA ) 

600 CONTINUE 
C 

IF ( KCM. EQ.O ) GO TO 300 
C 

IF (KSYMMZ. EQ.O ) NZL00P=2 
I F ( K S Y.MM Z • N E • 0 ) N Z LOOP = 1 
C 

NYLOOP=NY ( 1 ) 

C 

NX 1=NX ( 1 ) 

NY1=NY ( 1 ) 

NXY1=NX1*NY1 
DO 200 METH0D=1 , 2 
DO 170 I Y= 1 , NYLOOP 
DO 170 I Z= 1 , NZLQOP 
ITEMO=NXl*IY 
ITEM1=I TEMO— i 
ITEM2-1TEM0-2 
K Y= I Y +N Y 1 * ( I Z— 1 ) 

X1=PC ( 1, I TEMO) 
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153 

154 


C 


C 

170 

C 


190 


C 

C 

C 



199 


X2=PC ( 1 » ITEM1 ) 

X3=PCtl, ITEM2 ) 

XTE=0*5* (YPP( 1 » I T EMO ) +Y PM ( 1 » ITEMQ ) ) 

DET=(Xl-X2)*(X2**2-X3**2)-( X2-X3J * (Xl**2-X2**2) 

GO TO ( 153,152) , METHOD 

PHI 1*VELP0T( ITEMO J-+QEXP (UNI MAG*RED FRET BETA*PC( 1 , ITEMO? ) 

PH 1 2= VEL POT ( I T EMI ) *CE XP ( UN I M AG*RED FR E*B ETA4PC ( 1 » I T EMI ) i 
PHI3-=VELP0T( ITEM2)*CEXP(UNIMAG*REDFRE*BETA*PC( 1 , IT EM 2) ) 

GO TO 154 

PHI 1=VEL POT ( ITEMO i 

PHI 2=VELP0T ( I TEM1 ) 

PHI 3=VEL POT ( ITEM2) 

CONTINUE 

C0EF1=(+ ( PHI 1-PHI 2)*(X2**2-X3**2J- ( PHI 2-PH I 3 ) * ( Xl**2-X2**2 ) ) / 32 - r 
C0EF2-(- ( PHI 1-PHI 2) * ( X2-X3 ) + ( PHI 2— PHI 3 )*( X1-X2 ) ) /OET 
C0EF3-PHI1-C0EF1*XI-C0EF2’!‘X1**2 
VELPTE(KY)=OOEF3+COEF1*XTE+COEF2*XTE**2 
WRITE (NWRITE » 181 ) -PH 13, PH 1 2? PHI 1 , VELPTEl KY) 

I F ( METHOO • EQ • 2 ) 

IVELPTEIKY )= VEL PTE ( KY)*CE XP ( - UNIMAG -'‘RED FR ET8ETATXTE ) 

WRITE (NWRITE, 181) VELPTE IKY) 

CONTINUE 


CL=0. 

CM=0. 

AREAXY=0. 

ACH0RD=0. 

DO 190 I Y= 1 » NYLOOP 
ACL( I Y)=0 • 

ACM(IY) = 0. 

00 199 1 X= 1 , NX 1 

DO 199 I Y= 1, NYLOOP 

DO 199 I Z=1 , N ZLOOP 

1= I X+NX1* ( I Y- 1 ) +NXY 1* ( IZ-1) 

TWO EDGES OF EACH ELEMENT ARE ASSUMED TO BE PARALLEL TO X-AXIS 

DX=0.5*(YPP(1,I )+YPM( 1, I )-YMP( 1,1 ) -YMM ( 1 , I ) ) *BETA 

DY=0 • 5* ( Y P P ( 2, I )+YMP (2,1 >-YPM(2, I )-YMM(2,I ) > - 

DXY=DX*DY 

AREAXY=AREAXY+OXY 

I F ( KS YMMZ • EQ • - L ) COE FFT=4 • 

IF (KSYMMZ.EQ«+0)C0EFFT = 2.*{ 3.-.IZ*2) 

IF ( KSYMMZ • EQ. +1 )CQEFFT=0 . 


IF(KCM.EQ.1)AX1SMX=XLEZ+(XTEZ-XLEZ)*0.00 
I F ( KCM . EQ . 2 ) AX I SMX= XL EZ+ ( XT EZ-XLEZ ) *3 . 5 3 
I F ( KCM. EQ . 3 ) AX I SMX= 0. 5* ( CHAXI S ( I Y ) +CHAXI S( IY+1) ) 

IF (KCM. EQ. 4) AX I SMX-XLEZ-M XTEZ-XLEZ )*0.25 
ACL ( I Y ) =ACL( I Y ) +COEFFT* VELPGT ( I )*UN IMAG*REOFRE*DXY 
AC M ( I Y ) = ACM ( I Y ) +C0 E FF T * V EL POT ( I )*( l.-UNIMAG^REDFRE^T 8£TA*PC( 1 
l -AXISMX) )*OXY 

I F ( IX.NE.NXDGO TO 199 

ACHORD=ACHORD+ ( 0 « 5* (DCHORD ( I Y ) +D CHORD ( I Y+l ) ) )**2*0Y 
ACL ( IY)=ACLl IY) +COEFFT*V£LPTE( IY)*DY 

AC M ( I Y ) -ACM ( I Y ) -COEFFT*VELPTE( I Y)* ( 0. 5*( YPP ( 1 , I H-YPM(1 , 1 ) ! * 

1 B ETA-AX I SM X ) *DY 


CONTINUE 
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182 

181 

101 


102 

104 

195 

C 

107 


200 

300 

C 


c 

c 

10 

400 

J* * 

s> ■■ 


ACHQRD=ACHORD/AREAXY 

WRI TE (NWRITE » 182) AREAXY» ACHORD, AX I SMX» llZZi KCM 
FORMAT ( / //5X , 3F 9. 5 , 5X, 51 6 ) 

FORMAT (5X»2F 9. 5, 3X» 2F9. 5,3X , 2F9.5, 3X #2F9. 5.3X, 2F9. 5) 
WRITE(NWRITE.lOl) 

FORMAT ( / / ) 

DO 195 1 Y=1 fNYLGOP 

WRITE! N WRITE, 1 0 2 ) AC L ( I Y ) » AC M ( I Y ) 

ABSCL=CA BS ( ACL ( 1 Y ) ) 

ABSC M= CABS (ACM ( IY) ) 

REALCL=REAL { ACL (IY) ) 

REA.LCM-REAL ( ACM ( I Y) ) 

AIMGCL-A IMAGt ACL ( IY)) 

A IMGCM-AI MAG (ACM ( IY ) I 

FS AGCL= ATAN2 ( A I MGCL fREALCD^lSO./ 3.14159 
FSAGCM=AT A.M2( A I MGCM ,RE ALCM ) * 18 D. /3 . 14159 
WRI TE (NWRITE , 10 2 ) AB SCL t F SAGCL , ABSCM, FSAGCM 
F0RMAT(5X,2E15. 5,4X, 2E15.5) 

WRITE (NWR ITE » 104 ) 

FORMAT (/ ) 

CL=CL+ACL(IY)/AREAXY 
CM=CM+ACM( IY ) / ( AREAXY*ACHORD ) 

CONTINUE 


WRITE (NWR ITE, 107 )CLtCM»REDFRE 

F0RMAT(///5X. 'TOTAL LIFT COEFFICIENT - » 2E1J.5//5X 

1* TOTAL MOMENT COEFFICIENT = • ,2E15. 5. 5X, 'FREQUENCY 
ABSCL=CABS ( CL ) 

ABSCM=CABS ( CM ) 

REALCL=REAL ( CL ) 

REALCM=REAL(CM) 

A I M G C L = A I M A G ( C L ) 

A1MGCM=A IMAG (CM ) 

FSAGCL=ATAN2( AIMGCL .REALCL ) *180. /3. 14159 
FSAGCM= AT AN 2 ( AIMGCM»R EALCM )* 180. /3. 141 59 
WRITE (NWR ITE, 107 ) ABSCL , F SAGCL , ABSC M, FSAGCM , REDFRE 


» F7 .3// ) 


CONTINUE 

CONTINUE 


DO 400 I S= 1 , NS 
NXX=NX(IS) 

NYY-NY ( IS ) 

DO 10 IY*1,NYY 
DO 10 I X=1 , NXX 
I F ( 1X.EQ.NXX1G0 TO 10 
1 = l X+NX ( I S-)*( IY-D + NTI IS) 


J=I + 1 

DDX=PC( 1 , J )-PC( 1, I ) 

XPCM=0. 5*( PC ( 1 , J) +PC( 1, I ) ) 

SOURCE ( I ) = - 2 • 4 CEXP l-'JNI MAG^REDFRE^BET A^XPCM ) * 

( SOURCE ( J )*CEXP( SGNEXP*UNIMAG*REDFRE*PC( 1 , JJ/3ETA) 

1 -SOURCE ( I )#CEXP(SGNEXP*UNIMAG s ! t REDFRE*' 5 C(l , I ) / BETA) ) / 

Source ( i )=-2l*( ( velpoti j)-velpotit ) )/(beta*odx)+o.5*unimag*redfre 
l # (VELPOTI J ) +VELPQT( I )) ) 


CONTINUE 

CONTINUE 
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IF«KPRINTC8J .EQ.ll CALL PRI NTB (NTOTAL t 1 , AA, SOURCE, VELPOT , A ) 
I F ( KPRI NT (9) • EQ • 1 ) CALL PRINTB(NTOTAL , 1 ,AA, SOURCE* VELPOT , 5 ) 
I F ( KPRI NT ( 10 ) t EQ«. 1) CALL PRI NTS! NTOTAL* l » AA, SOURCE* VELPOT ,6 ) 

RETURN 

END 


SUBROUTINE. SOLUTNl NTOTAL , NT2S , AA, SOURCE t VELPOT ) 

COMMON/Z ZZ1/NX [ 20 ) , NY ( 20 ) ,NXY( 20) t NWt KSYMMY , KSYMMZ tNSYMMY , NSYMMZ 
COMMON/ ZZZ2/T AU , SPAN, TANGLE , TANGTE , CHORD, I ZZZ ,U 4ACH,REFLEN 
COMMON/Z ZZ99/ KPRI NT (10) , NREAD, NWRI TE 
DIMENSION AA ( NT2S )» SOURCE ( NTOTAL )» VELPOT ( NTOTAL ) 

IF(KPRINT!5).EQ.l)CALL PRINTS! NTOT AL, NT2S »AA , SOURCE , VELPOT , 1 ) 

IF (KPRI NT (6) . EQ.DCALL PRINTS! NTOTAL ,NT2S,AA, SOURCE, VELPOT, 2) 
T0L=J.00l 

CALL CCGELG! SOURCE, AA , NTOTAL , 1 ,TOL » I ER ) 


IF! I ER.NE.0)WRITE(N WRITE *101)1 ER 
IF! IER.NE.OJCALL DEBUG! 1001) 

FORMAT ( / // /2X , 1 — I ER = SK,' 

IF (KPRI NT (7) . EQ.DCALL PRINTB! NTOTAL, NT 2S, 
RETURN 


•///) 

AA, SOURCETVELPOT ,3 


) 


END 



SUBROUTINE PRINTB! MTOTAL , NT2 S, AA, SOURCE , VELPOT ,NPR I NT) 

COMPLEX SOURCE, VELPOT, AA, UNIMAG 

COMMON/ Z ZZ 1 /NX { 20 ) * NY ( 20 ) , NXY( 20) , NW , KSYMMY , KSYMMZ .NSYMMY , NS YMMZ 
C0MM0N/ZZZ2/T AUrSPAN, T ANGLE , TANGTE « CHORD , I ZZ Z t UMACH , RE FL EN 

CQMM0N/ZZZI1 / AL FA» ALFABC » REOFRE ♦ BOOYRf XLEZ * XTEZ» XNOSE. XTAI Lt KCM 

C0MM0N/ZZZ12/NS FX , NSBODY, NS , NT! 20 ) » KNORML ( 20 ) * KOI AF ( 20 ) * IS FACE ( 20 ) 
C0MM0N/ZZZ99/K PRINT (10) ,NREAD,NWRI TE 
C 

DIMENSION SOURCE! NT OTAL ), VELPOT (NTOTAL )» AA! NT2S) 

DIMENSION ABSVAL! 100) ,F ASEAN l 100) 

NY4=4* (NY! 1 )-l) 

GO TO! 1 » 2 , 3 * 4, 5, 6) , NPRINT 

1 CONTINUE 

WRITE INWRITE, 110) 

110 FORMAT ( ///2X, ' DISTRIBUTION OF AA(I ,J) '/) 

DO 11 1=1, NTOTAL 

WRITE (NWRI TE , 111)1 
N 1= I 

N 2 = NT OT A L * NT OTAL 

111 FORMAT C2X, ' INDEX=' , 12) 

WRITE (NWRI TE,1 12 )( AA! K) ,K=Ni,N2, NTOTAL) 

112 F0RMAT(1X,8E15,6) 

11 CONTINUE 

RETURN 

2 CONTINUE 

WRITE INWRITE, 221) i 

221 FORMAT ( ///2X , 'THE DISTRIBUTION OF SOURCE') 

225 CONTINUE 

NSBT0T=0 , • 

DO 229 I S= 1, NS 
I F ( NX Y ( I S ) • E Q • 0 ) GO TO 229 
WRITE (NWRI TE, 223) I SPACE! IS) 

223 FORMAT !//5X, • FOR PART *,I3) \ 

I NO=NSBT OT . 

NSBTOT=NSBTOT+NXY(IS) ! 

I FIN=NSBTOT ' I 

NXX=NX( IS) ’ j" 

DO 226 I X=1 , NXX • • 

WRITE! N WRITE, 228) . j 

IFINPRINT.GE.4.AND. (IX.EQ.NXX) )G0 TO 226 | 

IND=IND+1 , . I 

WRITE (NWRI TE, 227) (SOURCE ( KK) ,KK=IND, IF IN, NXX) | 

226 CONTINUE E 

227 FORMAT ! IX , 8E15 • 5) 'f 

228 FORMAT!/) ' • l 

229 CONTINUE \ 

RETURN p 

3 CONTINUE' 

W R I TE ( N WR I T E , 3 30 ) L 

330 FORMAT (Z//2X, 'THE DISTRIBUTION OF THE VELOCITY POTENTIAL ) 

GO TO 225 

4 CONTINUE 

WRITE ( NWR I TE , 440 ) 

440 FORMAT ! ///2X » ' THE DISTRIBUTION OF CP’) i 

GO TO 225 - i 

5 CONTINUE j 

WRITE! N WRITE, 550) j- 

i 
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550 F0RMAT(///2X, 'DISTRIBUTION OF ABSOLUTE VALUES AND PHASE-ANGLES • 
1,'OF CP' ) 

DO 555 KK=1» NTOTAL 
ABSVAL(KK)=CABS(SOURCE<KK) ) 

REA LCP= R E AL ( SOURCE ( KK ) ) 

AI MGCP=A I MAG ( SOURCE (KK )) 

FASEANtKK )=90. 

IF ( REALCP.NE.O « ) FASEANt KK )=ATAN2( A I MGCP * REALCP )*189./3. 14159 

555 CONTINUE 

551 CONTINUE 
NS3T0T=0 

DO 559 I S= 1 » N S 

IFINXYt IS) .EQ.0)G0 TO 559 

WRITE (NWRITEt 553) ISFACEI IS) 

553 FORMAT { / /5X » ' FOR PART • .12) 

I ND=NSBT OT 

NS8T 0T-=NS3T0T+NXY (IS) 

I F IN=NSBTOT 
NXX=NX( IS) 

DO 556 IX=1,NXX 
WRITE (N WRITE t 558) 

I F( NPRI NT.GE. 4. AND. ( I X. EQ.NXX) ) GO TO 556 
IND=IND+1 

WRITE(NWRITE,557) (ABSVAL(KK) , KK= I NO , I F IN ,NXX ) 

WRITEtNWRITE, 557) (FASEANt KK) , KK=I ND , I F IN , NXX ) 

556 CONTINUE 

557 FORMAT (1X»8E15.5) 

553 FORMAT (/) 

559 CONTINUE 
RETURN 

6 CONTINUE 

WR I TE ( NWRITE » 60 1 ) 

601 FORMAT ( ///2X» 'THE DISTRIBUTION OF CORRECTED VELOCITY POTENTIAL'/) 
DO 600 1=1, NTOTAL 
600 SOURCE! I ) =VEL POT ( I ) 

GO TO 225 
ENO 
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c 

c 

SUBROUTINE CCGELG(R,A,M,N,EPS, IER) 

C 

C 

COMPLEX WAKE, SO, R» PIVI ,TB,A 
DIMENSION A( 1 ) ,R( l) 

I F ( M ) 23, 23,1 
C 

C SEARCH FOR GREATEST ELEMENT IN MATRIX A 

1 IER=0 
PIV=0. 

MM=M*M 

NM=N*M 

DO 3 L=1 , MM 
TBB=CABS ( AID ) 

I F ( TB B -P I V ) 3 » 3 » 2 

2 P I V=T 8B 
I = L 

3 CONTINUE 
T0L=EPS*P1V 

C A(I) IS PIVOT ELEMENT. P1V CONTAINS THE ABSOLUTE VALUE OF A(I). 

C 

C 

C START ELIMINATION LOOP 

LST = 1 

DO 17 K-lj M 
C 

C TEST ON SINGULARITY 

IF(PIV)23,23,4 

4 I F ( IER)7»5,7 

5 IF(PIV-T0L)6,6,7 

6 I ER=K— 1 

7 PIVI=1./A( I ) 

J=( I-D/M 
I=I-J*M-K 
J=J+1-K 

C 1+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT 

C 

C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R 

DO 8 L-K , NM , M 
LL=L+ I 

TB=PIVI*R(LL) 

R ( L L ) =R CL) 

8 R(L)=TB 
C 

C IS ELIMINATION TERMINATED . 

IF(K-M)9, 13,18 
C 

C COLUMN INTERCHANGE IN MATRIX A 

9 LEND=LST +M-K 
IF( J) 12* 12, 10 

10 I I = J-^M . 

DO 11 L=LST , LEND 

T B-A( L ) 

LL=L+I I 
A I L ) =A( LL ) 

11 A ( LL ) =TB 
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C ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A 

12 DO 13 L=LST » MM* M 
LL=L+ I 

T B=P I VI #A ( LL ) 

A { LL) =A ( L ) 

13 A(L)=TB 
C 

C ' SAVE COLUMN INTERCHANGE INFORMATION 
A( LST ) =J 

C • 

C ELEMENT REDUCTION AND NEXT PIVOT SEARCH 

PIV=0. 

LST=LST+1 

J=0 

DO 16 I I =LST , LEND 
PIVI=-A( II ) 

I ST = 1 1 +M 
J = J + 1 

DO 15 L= I ST * MMf M 
LL=L-J 

A(L)=A(L)+PIVI*A(LL) 

TBB=C ABS (AIL) ) 

1FITBB-P IV) 15, 15, 14 

14 P IV=TB8 
I =L 

15 CONTINUE 

DO 16 L=K,NM,M 
LL=L+J 

16 R(LU=R(LL )+PIVI*R(L) 

17 LST=LST+M 
PND OF ELIMINATION LOOP 


BACK SUBSTITUTION AND BACK INTERCHANGE 

18 I F ( M- 1 ) 23 , 22 * 19 

19 IST=MM+M 
LST =M+1 
DO 21 1=2, M 
Il-LST-I 
I ST-1ST— LST 
L-IST-M 
L= A C t. 

00 21 J= I I , NM , M 
TB=R(J) 

LL=J 

DO 20 K= I ST , MM , M 
L L = L L + 1 

20 TB=TB-A { K ) *R ( LL ) 

K= J+L 
R( J)=R(K) 

21 R ( K ) =TB 

22 RETURN 
C 

. C 

c ERROR RETURN 

23 I ER=- 1 
RETURN 
END 
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C 


SUBROUTINE COEFF ( NTOTAL » PC » YPP » YPM . YMP * YMM *KWAKE t NT2S» A A , SOURCE, 
ISCTRAN ) 


THIS SUBPROGRAM IS TO CALCULATE DOUBLET AND SOURCE FOR SUBSONIC 

FLOW 


COMPLEX SOURCE ♦SCTR AN tAA, UNIMAG, XNUNST tCMPXDB f CMPXSC,UIOMER UNST 

COMPLEX XNTRAN 


C 


C 

c 


COMMON/ 2 2Z1 /NX { 20 ) » NY ( 20 ) » NXY ( 20 ) » NW» KSYMMY , K SYMMZ »NSYMMY * NSYMMZ 
COMMON/ Z2Z2/TAU , SPAN, TANGLE , TANGTE , CHORD » I ZZZ, UMACH, REEL EN 
COMMON/ZZZIO/K WAKES (20 ) , NWAKE ,KXI NCR 

COMMON/Z / Z1 1/ ALFA* ALFABC » REDFREf BODYR « XLEZ * XTEZ * XNOSEi XTAIL.KCM 
COMMON/ ZZ Z 12/ NS FX , N S80DY « NS * NT ( 20 ) * KNORML ( 20 ) , KD I AF ( 20 J , I S FACE ( 20 ) 
COMMON/ ZZZ13 /AC HORD ( 10 ) »XCHORD( 10) 

COMMON/ ZZZ88/QV ( 3 ) » QQ t I . J » I SYMMY, I SYMMZ 
COMMON/ Z ZZ99/KPRI NT ( 10) , NREAD , NWR I TE 

DIMENSION YPP ( 3 » NTOTAL ), YPM ( 3 , NTOTAL ) , YMP < 3 , NTOTAL ) , YMM ( 3 , NTOT AL ) , 
1PC (3 iNTOTAL) , SOURCE (NTOTAL ) , SCTRA'N ( NTOTAL ) ,K WAKE ( NTOTAL ) , A A ( NT2S ) 
DIMENSION A1 ( 3 , 2) , A 2 ( 3 , 2 ) * A 1A1 ( 2 ) * A2A21 2) , QCRA 1 ( 3 ) » QCR A2 ( 3 ) 
DIMENSION A1V(3) , A2V( 3) , YNORM (3 ) , UN ( 3 ) , SN ( 3 ) , S NU N ( 3 ) „ A 1 A 2 ( 2 » 2 ) 
DIMENSION PZ(3) »Q ( 3 * A ) » OMCR A 1 ( 3, 2 ) » QMCRA21 3 » 2 ) » A A 1 { 3 , 2 ) * A A2 ( 3 » 2 ) 

D I ME NS I ON VMM (3 ) » VDD( 3 ) t VMD ( 3 ) » YTEPZt 3 ) ♦ YTEMZl 3 ) » SIGNPT ( 3) 

DIMENSION AVA1I3) » AVA21 3 ) » A1CRA2( 4 ) * QM ( 3 ,4) 

DIMENSION WPP 13) »WPM( 3) ,WMP(3) , WMM l 3 ) , WPC(3) 


D0TPR0(X1,Y1»Z1,X2,Y2,Z2) = Xl*X2-*-Yl*Y2+Z 1^Z2 

PROMIX ( XXI, YYl,ZZl,XX2,YY2,ZZ2,XX3,YY3,ZZ3)*(YY2*ZZ3~YY3*ZZ2)*XXi 
i -(XX2*ZZ3-XX3*ZZ2)*YYU-(XX2*YY3-XX3*YY2)*ZZ1 


C 


2 

C 


9911 

9922 

9933 


IF (UMACH. GT.l.lCALL DEBUG (400) 

BETA=SQRT( l.-UMACH**2) 

WRITE(NWRITE,2)NSB0DY,NS, (NT(I ) ,1 =1 ,NS) » ( KNORML (J),J-1»NS)»(KDIAF 
1 ( K ) , K= 1 , NS ) »( I SPACE (I I ),I 1=1, NS) , ( NXY ( KK ) , KK= 1 . NS) ,NSYMMY , NSYMMZ 
F0RMAT(3X,25I5) 


IF ( IZZZ.GE. 100. AND* I ZZZ .LT. 200) WRI TE (N WRITE »9911 ) 
I F C IZZZ.GE. 200 .AND. IZZZ.LT. 300) WRITE ( N WRITE » 9922 ) 
I F ( IZZZ.GE. 300. AND. I ZZZ .LT.400 ) WRI TE (N WR ITE » 9933 ) 


FORMAT l // 2X , • OSCILLATION IN FIRST BENDING MODE 

FORMAT ( / /2X , OSCILLATION IN TRANSLATION 

FORMAT ( / / 2X » • OSCILLATION IN PITCH ' 


IF CI2ZZ.EQ.301) XP ITCH=XLEZ + ( XTE Z-XLE Z ) *0 .00 
IFUZZZ.E0.302) XP ITCH=XLEZ+ ( XTEZ-XLEZ) *0. 50 
IF(KCM*EQ. 11 ) X P lTCH=XLEZ+( XTEZ-XLEZ J*0. 15 
I F ( KCM.EQ. 12)XP I TCH = XL EZ+ (XTEZ-XLEZ ) *0 . 20 
IFtKCM.EQ. 13) XP I TCH=XLEZ+( XTEZ-XLEZ) *0.25 
IF (KCM.EQ. 14 )XPITCH=XLEZM XTEZ-XLEZ) *0.30 
IF (KCM.EQ. 15 >XPITCH=XLEZ+( XTEZ-XLEZ )*0. 325 
IF ( KCM. EQ. 16) XP I TCH=XLEZ+( XTEZ-XLEZ )*0. 350 
1F( KCM.EQ. 17) XP I TOH=XLEZ+( XTEZ-XLEZ ) *0. 1875 
IF( KCM.EQ. 13JXP I TCH=XLEZ+( XTEZ-XLE Z ) *0 • 4500 



UNST 


UNSI 


UNST 

» //) 

UNST 

«//) 

UNST 

'//) 

UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


UNST 


C 
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c 


c 


UNI MAG= ( 0 . » 1 . ) ™ 

OMEGA=REDFRE*UMACH/BETA UNb 

ALFRBC=ALFABC*3. 14159/180. 

S I N ABC=S I N l ALFRBC ) 

COS ABC=C0S ( ALFRBC ) 

00 6 NN=1*NT2S 
A A { NN ) = 0 • 

INSERT THE KRONECKER DELTA ISWITCH FOR THE LOWER SIDES OF DIAPHRAGMS) 


DO 9 IS= 1 1 NS 
NXY S=NXY (IS) 

DO 9 I 1 = 1 ? NXY S 
1 * 1 I+NTl IS ) 

J = I 

NNM=I+( J-l ) ^NTOTAL 
AA ( NNN ) = 1 . 

CONTINUE 


DO 10 1=1 t MTOTAL 
SOURCE 1 1 ) = 0 • 
SCTRAN ( I ) =0 . 

10 CONTINUE 
C 

SIGNPT ( 1 ) = 1.0 


C 

CONST=+C .5/3. 14159 
C 

DO 180 JFACE=lf NS 
J SBTQT=NXY ( J FACE ) 
DO 180 J J = 1 » J SBTOT 
J=NT( JFACE)+JJ 


C 


1002 


1003 

C 


100 4 
C 


DO 1002 K= 1 » 3 

AA1(K, 1)=0. 5*(YPP(K, J )-YMP(K» J ) ) 

AA1 ( K » 2 ) = 0 • 5* ( YPM ( K » J ) - YMM( K * J ) ) 

AA2(K, 1)=0.5*<YPP(K,J )-YPM(K f J ) ) 

AA2(K,2)=0. 5*<YMP(K, J )-YMM(K» J)) 

A 1 ( K » 1 ) = A A 1 1 K ? 1 ) 

A1 1 K» 2) =AA1 1 K ? 2 ) 

A2(K» 1) = AA2(K, 1) 

A 2 ( K? 2 }=AA2 1 K » 2) 

CONTINUE 

A1A1 ( L)=DOTPRO l A1 1 l.?L ) * A1 ( 2»L ) » A 1 ( 3 » L ) * All 1? L) ? A1 ( 2tL) » A1 13 » L ) ) 

A2A2 ( L) = DOTPRO( A2( 1 «L ) » A2 ( 2 » L ) » A2( 3* L) » A2( 1»L) » A2 l 2* L) t A2( 3 r L ) ) 


DO 1004 K= 1 * 3 

IF( A1AK 1 ) .EQ.O.JAKK, 1)=A1(K» 2) 
I F t A 1 A 1 ( 2 ) • E 0 • 0 • ) A 1 ( K » 2 ) = A 1 ( K » 1 ) 
IF ( A2A2 ( II • EO »0 • ) A2 ( K » 1)=A2(K?2) 

IF(A2A2(2).E0.0.)A2(K,2)=A2(K,1) 

CONTINUE 


DO 1006 L= 1* 2 

A1A1C L) = DOTPRO( All ItL > * Al(2jL) t Alt 3»L) > All I* LI ? Alt 2? L) ■, All 3» L ) ) 


SUB 
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1006 A2A 2 ( L) = DOT PRO ( A2( l»L>»A2(2*L)»A2(3,L)»A2(l*L)fA2(2»L)»A2(3,L)) 

C 

DO 1007 L= 1» 2 
DO 1007 M=1 * 2 

1007 A1A2(L»M >=DOTPRO( A1 ( 1 f L ) » All 2 »L ) » AH 3 » L ) » A2( 1 » M) ♦ A2 ( 2* M) « A2( 3* M) I 
C 

A1CRA2I 1)=SQRT(A1A1(2)*A2A2( 1)-A1A21 2* 1)**2) 

AICRA2 { 2 )= SORT ( A1A1 ( 1 )^A2A2( 1 ) -A1A2( 1* 1)**2) 

A1CRA2(3)=SQRT( A1A1( U*A2A2<2)-A1A2( l v 2)**2) 

A1CRA2J 4>=SQRT( A1A1< 2) *A2A2 (2J-A1A2 ( 2 .21**2) 

C . 

DO 1008 K=l,3 

AVA11K) = 0.5 1 MA1(K»1KA1(K>2)) 

1008 AVA2(K)=0.5*< A2 ( K ♦ 1 )+A2 ( K ♦ 2 ) ) 

C 

YN0RM(l)=AVAl(2)*AVA2(3)-AVAi(3)*AVA2(2) 

YNQRMI2)— AVA1I2 )*AVA2(1 )-AVAl ( 1 )*AVA2( 3 ) 
YNORM(3)=AVAi(l)«AVA2(2)-AVAl(2)*AVA2( 1) 

C 

SN( 1|=Y NORM (1)/8ETA 
SN( 2 )=YNORM ( 2 ) 

SN ( 3 )=YN0RM ( 3) 

ASN=SQRT ( SN { 1 )**2.+SN( 2)**2+SN(3)**2) 

AYN=SQRT (YN0RM( 1)**2+YNQRM( 2)**2+YNORM(3)**2) 

DO 1010 1 * 3 

IJN(K) = YNORM(K)/AYN 
SNUN(K)=SN(K)/ASN 
1010 CONTINUE 
C 

c . 

DO 170 I FACE= 1 » NS 
C 

I SBTOT=NXY< IFACE) 

DO 160 1 1 = I » ISBTOT 
I =NT ( I FACE ) + I I 
C 

00 160 I SYMMY= 1* NSYMMY 
DO 160 ISYMMZ=1 ,NSYMMZ 
SI GNPT ( 2 ) =3 •— 2* I SYMMY 
SIGNPT(3)=3.-2*ISYMMZ 
C 

DO 1102 K=1.3 

1102 PZ (K >=PC ( Kt J )-PC ( K, !I*S1GNPTIK) 

QOOTUN=DOTPRO(UN( 1) ,'JN(2) ,UN(3),PZ(1) , PZ(2) , PZI3 ) ) 

C 

DO 1110 K- 1 * 3 

0 ( K i 1 ) =Y PM ( K » J )-PC ( K» I) *S I GNPT ( K ) 

Q(K,2)=YPP(K,J)-PC(K, I )*SIGNPT(K) 

Q(K,3)=YMP(K, J)-PC(Kt I )*SIGNPT(K) 

Q ( K * 4 ) =YMM( K t J ) -PC ( Kt I ) *S I GNPT (K ) 

1110 CONTINUE 

DO 1111 K=1 1 3 

QM ( K 1 1 ) =0 . 5* ( 0 ( K t l ) +Q 1 K » 2 ) > 

QM(K,2)=0.5*(Q(K,2)+0(K.3) ) 

QM(K,3> = 0.5*(Q<K,3)+Q(K,4> > 

QM ( K , 4 ) =0 . 5* ( Q ( K , 4 ) +0 ( K » 1 ) ) 

1111 CONTINUE 
C 
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1112 

C 

C 

5502 

5504 

5506 

5508 

5510 

C 


5520 

C 

C 

C 

C 


DO 1112 K= 1*3 
KP1=K+1 

KP2=K+2 

l F ( KP 1 • GT • 3 ) K P1=K# 1-3 
IF(KP2.GT.3)KP2=KP2-3 

QMCRAKK, IJsQM KPl » 2) *A 1 ( KP2 , 1 )-QM ( KP2 , 2) *A 1 (KPl » 1 ) 

0 M C R A 1 < K , 2 ) = Q M ( K P 1 , 4 ) * A 1 ( K P 2 , 2 ) - Q M ( K P 2 , 4 ) * A 1 ( K P 1 , 2 ) 

QMCRA2(K,1)=QM(KP1, 1 ) *A2 ( KP2. 1 )-0M ( KP2 , 1 ) *A2 ( KPl , 1 ) 

QMCRA2(K , 2 ) =OM ( KPl , 3) *A2 ( KP2 , 2) -QM ( K ?2 , 3 ) *A2 ( KP 1 » 2 ) 

CONTINUE 

RC=SQRT (DOTPRO(P2( 1 ) , PZ ( 2 ) , PZ ( 3) * PZ( 1) » PZ( 2) . PZ( 3) ) ) 

U I OMEP.-UNI MAG*OMEGA*RC 
CMPXDB=CEXP ( - U T OM ER ) *( i.+UIOMER) 

CMPXSC=CEXP(-UIOMER)*CEXP(-UNIMAG*QM£GA*UMACH*PC(l ,J) ) 

DO 155 IC0RNR=1,4 

GO TO (5502,5504, 5506,5508) , ICORNR 
CONTINUE 
$ I GN 12=-1 • 

ICSI=1 
I ET A=2 
GO TO 5510 
CONTINUE 
SIGN12=+1. 

ICSI=1 
I ETA = 1 
GO TO 5510 
CONTINUE 
S1GN12=-1. 

I CSI=2 
I ETA=1 
GO TO 5510 
CONTINUE 
SIGN 12-+ 1 • 

I C S I = 2 

I ETA=2 • "*> 

CONTINUE 

DO 5520 K= 1,3 
QV(K)=Q(K, ICORNR) 

A1V ( K ) =A 1 ( K, IETA ) 

A2V ( K) =A2 ( K, ICSI ) 

QCRA1 (K)=0MCRA1(K, IETA) 

QCRA2 ( K ) =QMCR A2 (K , ICSI ) 

CONTINUE 

QQ=SQRT( OOTPROt GV( 1 ) , QV { 2) , QV { 3 ) , QV ( 1) ,QV( 2) »QV{ 3)1 ) 

CALL LOG ( I CORNR , A IV ,QCRA l , AL0G1 , 1 ) 

CALL LOG ( I CORNR , A2V ,QCR A2 , AL0G2, 2 ) 

T A N P - 0 . 

I F ( QDOTUN. EQ . 0 • ) GO TO 5550 

BNUMER=-D0TPR0( OCRAK 1) » OCR A 1 ( 2 ) , OCRAl ( 3 ) ,QCRA2( l) ,QCRA2(2>, 

1 OCR A2 ( 3 ) ) 

DEN0M=QQ*3D0TUN*A1CRA2{ ICORNR) 


UNST 

UNST 

UNSI 

UNST 


SUB 


SUB 
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IF C DENOM.NE.O. )TANP=ATANP ( HNUMER, DENOM) 

5550 CONTINUE 
C 

COEFFI=DQTPRO(UN(l) ,UN( 21 ,UN(3) ,QCRA 1 ( 1 ) , QCRA1 (2 ) ,0CRA1( 3) ) SUB 

COEFF2=DOTPRO(UN( 1 ) ,UN ( 2) ,UN ( 3) ,QCRA2( 1 ) , QCRA2 ( 2 ) , QCRA2 ( 3) ) SUB 

C 

SRCINT=+CONST*SIGN12* SUB 

1 (-COEFFl*ALOGL+COEFF2*ALOG2-QDOTUN*TANP) 

C 

DBTINT=-CONST*SIGN12*TANP 

C 

SSNINT*l. 

I F ( ISYMMY . EQ • 2 ) SGNI NT=SGN INT*KSYMMY 

I F ( I SYMMZ . EQ. 2) SGNI NT= SGNI NT*K SYMMZ 
C 

NNN=I ♦( J-l )*NTOTAL 
C 


AA ( NNN ) =AA ( NNN )— SGN INT*DBTINT*CMPXDB UNST 

C 

ARG=ABS ( (PC(2»J)-BODYR)/( 0. 5* SPAN- BOD YR ) ) 

I F ( I III . EQ . 10 1 > XNUN ST =- UN I MAG* SNU.NI (3 )* UNST 

1 O.18043*ARG+1.70255*ARG**2-1.13688*ARG**3+0.25387*ARG**4) UNST 

IF (1 ZZZ.GE.200.AND.IZZZ.LT. 300) XNUNST=+UNI MAG*REDFRE*SNUN ( 3 ) UNST 

l F{ I ZZZ.GE.300.AND.IZZZ.LT. 400) UNST 

1XNUNST=+ ( UNI MAG*( PC ( l , J ) *8ETA-XPITCH ) *REDFR E+l • 0 ) *SNUN (3 ) UNST 


C 

C FOLLOWING IS FOR FLUTTER 

IF( IZZZ.LT . 5 0 0 ) GO TO 499 

XNUNST=+ (UNI MAG*( PC ( 1 , J )*BETA-XPITCH )*REDFRE+i.O )*SNUN ( 3 ) 
XNTRAN=+UNIMAG*REDFRE*SNUN< 3) 

SCTRANt I )=SCTRAN( I ) +SGN I NT*SRC INT*XNTRAN*CMPXSC 
499 CONTINUE 
C 

SOURCEt I )=S0URCE( I ) +SGNINT*SRCINT*XNUNST*CMPXSC UNST 

C 

150 CONTINUE 

I F ( I .GT . 0 ) GO TO 155 
I F ( J . GT • 5 ) GO TO 155 

WRITE(NWRITE, 660) l , J, ISYMMY, I SYMMZ ♦ 1C0RNR , DBTI NT , SRC I NT » COEFFI , 
1AL0G1»C0EFF2. AL0G2, SOURCE ( I ) 

660 FORMAT ( IX ♦ 5 14 , 8E13. 4) 

155 CONTINUE 
160 CONTINUE 
170 CONTINUE 
180 CONTINUE 
C 

C CALL TIME (10) 

C 

C FOLLOWING IS THE EVALUATION OF WAKE COEFFICIENTS 

C 

FACTO R= 1.0 
IXYWAK^O 
C 

IFtKXINCR.FQ.20) X I NCR2=2 .0*CH0RD/N WAKE ; 

IF(KXINCR.EQ.21)XINCRZ=2.5*CH0RD/NWAKE 
I F ( K X I N C R * E Q . 2 2 ) X I N C R Z = 3 . 0 * C H 0 R D / N W A K E 
I F( KX INCR.EQ . 23 ) X INCRZ=4.0*CH0R0/NWAKE 
IFTKXINCR.EQ. 10 l ) XI NCRZ=0 . 04*CH0RD/FAC TOR 
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c 

188 

C 

C 

C 


C 

770 

C 

C 

c 

6601 

C 

C 

C 


1 F { KX I NCR • EG* 102 ) X I NCRZ=0 . 06*CH0RD/FACT0R 
1F(KXINCR.EQ. 103) XI NCRZ=0.08*CH0RD /FACTOR 
I FIKXINCR.EO. 104) XI NCRZ=0. 10*CHORD/FACTOR 

WRITE! 6, 188 I KXINCR » NWAKE» XINCRZ 
FORMAT C//5X., 217 » F14. 4/'/) 

XINCR=XINCRZ 

00 999 I XWAKE= 1 »NWAK£ 

X I NCR=X I NCR* FAC TOR 
XWAKE 1=XWAKE2 
XWAKE2=XWAK£1+XINCR 

00 888 J=itNTOTAL 
I F ( KWAKEt J ) . EQ .0 ) GO TO 888 

WPP(1)=YPP(1»J J+XWAKE2 
WPP(2)=YPP(2» J ) 

WPP(3)=YPP(3, J) 

WPM ( 1 ) =YPM ( L.J ) +X WAKE2 
WPM ( 2 ) =YPM ( 2 » J ) 

WPM(3)=YPM<3, J ) 

WMP(1)=YPP(1»J) +XWAKE 1 
WMPt 2) =YPP (2 » J ) 

WMP ( 3 ) =YPP ( 3 » J ) 

WMM( 1 )=YPM( 1 1 J 5 +XWAKE1 
WMM ( 2 ) =YPM (2 ♦ J ) 

WMM ( 3) = Y P-M ( 3 ♦ J ) 

DO 770 K = 1 » 3 

WPC (K )= ( WPP ( K ) +WPM ( K ) +WMP ( K ) +WMMI K ) ) /4. 

A1V ( K )=0 . 5* ( WPP ( K ) —WMP ( K ) ) 

A2V ( K ) =0 . 5* ( WM P ( K ) -WMM ( K ) ) 

CONTINUE 

AIVA1V=D0TPR0{ A IV ( 1 ) t A1V ( 2) » A IV l 3 ) T A1V ( 1 ) * A1V ( 2 ) * A IV ! 3 ) ) 
A2VA2V=D0T PRO ( A2V(1)»A2V(2) ,A2V(3) , A2V( 1 ) ,A2V( 2) ,A2V(3) ) 
A1VA2V=D0TPR0( A IV ( 1 ) * A1V (2) » A1V (3 ) » A2V ( 1 ) » A2V (2) »A2V (3 ) ) 

A1CSA2=SQRT{A1VA1V*A2VA2V-A1VA2V**2) 

Y.NORM ( 1 ) =A1V( 2) *A2V (3 )-A IV ( 3 ) *A2V ( 2 ) 

YM0RM(2) = A1V{3)*A2V (1 )-AlV( 1)*A2V< 3) 

YNORM (3)— A1V ( 1 )*A2V (2 )-AlV( 2 )*A2V( 1 ) 

ABN0RM=S0RT(YN0RM(1 )**2 +YNORM < 2 ) **2+YN0RM { 3 ) **2 ) 

DO 6601 K= 1*3 

UN ( K )■= YNORM ( K ) / ABNORM 

co'ntinue 

00 777 I=1,NT0TAL 
NNN=I +( J-l) *N TOTAL 

DO 777 I SYMMY=L,NSYMWY 
DO 777 I SYMMZ = 1 » NSYMMZ 
S IGNPTl 2 ) =3.-2* ISYMMY 
S I GN P T ( 3 ) = 3 . - 2* I SYMMZ 

00 6602 K=l»3 
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6602 PZ(K)=WPC(K)-PC(K, I1*SIGNPT(K) 

Q00TUN-0QTPR0 ( UN ( 1 ) , UN ( 2 ) , UN{ 3 ) , PZ ( 1) , P Z ( 2 ) , PZ ( 3 ) ) 

C 

DO 6610 K= 1,3 

Q ( K , 1 ) =WPM ( K ) -PC ( K , I ) * S 1GNPT ( K ) 

Q(K,2)=WPP(K)-PC(K, I )*SIGNPT ( K ) 

Q(K,3J=WMP(K)-PC(K, I J*SIGNPT ( K ) 
q ( « , 4) =WMM (K)-PC(K»I ) * S I GMPT t K J 

6610 CONTINUE 

00 6611 K = 1 » 3 

OM(K, 1)=0.5*(Q(K, l)+3(K f 2) ) 

QM ( K , 2 ) =0 • 5* ( Q ( K , 2 ) +0 ( K » 3 ) ) 

QM(K,3)=0.5*(Q(K,3)+Q(K,4) ) 

QM ( K » 4 ) = 0 • 5* ( 0 1 K » 4 ) +3 { K » 1 ) ) 

6611 CONTINUE 
C 

DO 6612 K-1,3 

KPl=K+l 

KP2=K+2 

IF (KP1.GT.3 )!<Pl=KPl-3 
I F ( KP 2 « GT . 3 ) KP 2 = KP2— 3 

QMCRAKK, 1)=JMC KP1, 2) *A IV (KP2 ) -QM ( KP 2 , 2 ) *A IV ( KP 1 ) 

QMCRAUK, 2 ) = QM ( KP 1 , 4 > * A1 V ( KP2 ) -QM ( KP 2 » 4 ) *A1 V< KP1 ) 

QMCRA2 ( K » 1 ) =QM ( KP 1 » 1) *A2V { KP2 ) -QM { KP2 , 1 ) *A2V ( KP 1 > 

0MCRA2(K,2)=0M( KP1, 3 ) *A2V ( KP2 J -QM { KP2»3 )*A2V(KP1) 

6612 CONTINUE 
C 

RC = SQRT ( DOTPRO ( PZ ( 1 ) » PZ I 2) » PZ ( 3 ) , PZ( 1 ) , PZl 2 ) ,PZ(3) ) ) UN 

WAKE=0. 

C 

00 666 I C0RNR = 1 ,4 

GO TO (7702,7704,7706,7708) ,ICORNR 
7702 CONTINUE 

SIGN 12=— 1 • 

ICSI=1 
I ET A= 2 
GO TO 7710 
7704 CONTINUE 

SIGN12=+1. 

rcsi=i 

I ETA= 1 
GO TO 7710 
7706 CONTINUE 

S I GN 1 2 =- 1 . 

ICS 1 = 2 
I ETA= 1 
GO TO 7710 
7708 CONTINUE 

SIGN12=+1. 

I CS I =2 
I ETA=2 

7710 CONTINUE 
C 

DO 7720 K= 1 » 3 
QV(K)=0(K, ICORMR) 

QCRAKK) =QMCRA1(K, IETA) 

Q C R A 2 ( K } = Q MC R. A 2 1 K , ICS I ) 

7720 CONTINUE 
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SUB 


C 

C 


7750 

C 

c 

666 

C 


2244 

C 

C 


C 

C 


777 

888 

999 


QQ=SQRT ( D0TPR0 ( QV ( 1 ) » QV ( 2 ) » QV ( 3) * 0V 1 1) «QV( 2) «QV(3 ) ) ) 

TANP=0. 

IFIQDOTUN.EQ.0. )G0 TO 7750 

HNUMEP=- DOT PRO ( QCRA 1 ( 1 ) , OCR A 1 ( Z ) , QCRA 1 { 3 ) , QCRA2 ( 1 ) , QCRA2 (2) , 
1 QCR A2 { 3 ) ) 

DENQM=QQ*QD0TUN*A1C SA2 

IFtOENOM.NE.O. ) TANP=ATANP { HNUMER , DENOM) 

CONTINUE 

WAKE=WAKE-SIGN12*TANP*C0NST 
CONTINUE 
I XYWAK=I XYWAK+1 

I F ( IXYWAK.EQ. 1) WRITE! 6,2244 )I , J, ISYMMY, I SYMMZ, WAKE 
FORMAT (41 5, El 6. 5) 


SGNINT=1. 

I F C ISYMMY. EQ. 2) SGNI NT=SGNINT*KSYMMY 
I Ft ISYMMZ.EQ.2) SGNI NT=SGN I NT*KSYMM Z 

XXX= WPC ( 1 )-PC( 1 » J ) 

AA ( NNN ) = AA ( NNN ) -W AK E*SGN I NT* ( 1 . +UNI MAG*OMEGA *RC J * 

1 C EXP ( -UN I MAG*REDFR E* ( XXX+UMACH*RC ) /BETA ) 

CONTINUE 
CONTINUE 
CONTINUE 
RETURN 
END 


SUE 
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o o 


c 

c 

SU3R0UT INE LOG{ ICORNR * AV» QCRAV » ALOG * ICOOR) 
CQMMON/ZZZ83/QV (3 ) » 09» I » J f I SYMMY» I SYMMZ 
DIMENSION AV I 3 ) tQCRAV 13) 

C 

DOT PRO ( Xl»YltZltX2f Y2 t Z2)=X1*X2+Y1*Y2+Z1*Z2 

c 

AVAV=DOTPRO( AVI 1 ), AVI 2) • AVI 3 ) » AVI 1 ) , AV I 2 ) , AVI 3 ) ) 
QQAV=DOTPRG( QVI 1 ) * QV I 2 ) » QV I 3 ) , AV I I } , AV I 2} , AV 1 3 ) ) 
QXA=SQRT I DOT PRO! QCRAV I 1 ) » QCRAV I 2 ) , QCRAV I 3 ) » 

1 QCRAVI 1) ,QCRAV(2) ,QCRAV(3) )» 

ALOG= AS I NH! QQAV/QXA ) / SORT I AVAV ) 

RETURN 

END 


C 

FUNCTION ASINH(X) 

A8X=ABS I X ) 

XSQ=X*X 

IFI ABX.LE. 0,001 )ASINH=X*( 1.-0.1666666*XSQ+. 075*XSQ*XSQ ) 
I F l A BX • GT • 3 • 3 0 1 ) AS I NH = I X/ABX ) *ALOG I ABX-t-SQRT I l.+XSQ)) 
RETURN 
END 


FUNCTION ATANP I HNUM ER , DENOM ) 
ADENOM= ABS ( DENOM ) 

AT ANP=ATAN2 I HNUMER » ADENOM ) 
IFI DENOM. LT.O. ) AT ANP=-ATANP 
RETURN 
END 
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u c o u 


c 

c 

c 

r 

r 


2 

C 


qq* * 
qq? .? 
90 33 


CD 

rn 


SU°^PUTIMF CHFFFt K T n T ALtPC» YPP »'Y°M»VVP f YMMf SUN03 »NT2$»AA* SOU 3 CF * 
ISCT P AN ) 

r u t c; <;ijnppMG p AM IS T U CALCULATF DOUBLET 
A\n r,rtJ n CR C C° UNSTEADY SUPER S^Nir FLCW 

r CMP L FX F C UP ‘ E , SC T p AM t A A , U N I w A C , X NU N S T t X N TP AN , C y P X DB , C M d v SC 

CCMfH EX Cr^GA, SIMHUt CHSHU 

rrMMMN/ZZ7l/MX(P0),NV(P0),NXY(20),NN,KSYMMYASYMM7 tN <;YMMY,NSY^7 

r c"«rN/zz|'o/'t Ji«FSt’ci sc En , rocn , «nov» ,x. ez.xti-z .xhoxi;, *t»ii. c« 

rCHMrN/Z7ZWNKCTF,K->BtKT(*),HOMU«I4.?».0<P»->' 

-rMtnN/ZJZ77/ 2 V(3) .^.KORFm ,sr,WU4t2! ,'iI..NV2{A. . 

CCYMCN/7.7 733/ I ,.J, ISYj’YY, 1 SY«M?. 

r CMMrN/ZiZqD/KPPI N T (' 0) v'm(q K'TnT a i ) , YM« ( 3 , N T UTM ) » 

D0 Tpr ' 0( X 1 1 v l » 7. L ?X2 » Y2 i Z2 ) =X 1*X2+Y L -Y2 + 7 l "7.2 

SUFfffIXX* »YY* »?Z? . XX 7.YY’*ZZ?.1*TX* *XX2-YYl*YY2-Z Z1*ZZ2 


T F ( U f * ACH . ED . 1 . ) U w AfH- SQR T ( 7 . ) _ 

I F ( 'JN AC.H . L T . I • ) C ' LL nF p UFI500) ~ ; 

BETA=SQRTUJYArH*^?-l. ) , K N np M l ( j ) , J=’ , NS ) , ( KO I AF ( K } , 

^Iiritir.^uV.uIiUs.I.NXvUK, .XX.UNS, .NXYMMY.NSYHMZ ' 
FORMAT ( 3X, 3515) 

T c ( i ? ; 7 . DC . 100 , AND . 1 7 7 Z . L T . 230 ) W p ! T E ( .6 » 90 U ) 

;fI!z^:GF. 2 00. M C.I7Z2.LT.300;^.^J6,^22 

VsrJuiTIci li FI»^ P.gNW, - 

pop m.AT ( / / 7 X * I SC i L LA I 1L.\ I' ' . 

cnPMMf//?X,» — OSCII.LATICN TP ANSLA T ON - 

F^«t!//2x: RSCUUX.CN !« -ITCH 

X p I TCH=XL f 7 + ( X T E7-XLF 7. ) *^ TI ' vi-7W C q m 

• tc| 1 7 7 7 F p 3-3 L ) X° I T C,U = XL F7 f ( X T ?:Z-XL>: 7 ) vD.03 

IM!777:Fi:T5FIXPlTCH = XLr»>(XTFZ-XieZ>*J.51 
IFlir.«.FO.H)XP'TCH = )(LFZtlX-FZ-XI -FMTO.1^ 

I F ( K r M • ED • 1 2 ) X n T T CH=X LE 7. +■ ( X E Z XL -. ) • 


UNS T 


-—«//) 

UNS T 

UNS T 

UNST 

UN'S- 

• it ) 

UNS T 

'// ) 

JNS T 


UNS T 

UNST 

UNS 1 

UNST 
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CD 

CD 

CO 

CD" 

CD 

CD 

C 


C 


c 

6 

C 

C 

c 


c 

c 

c 


9 

c 

c 


in 

c 

c 


c 


IF(K r M.EQ.l3)XPI' r CH=;XLEZ + (XTFZ— XLEZ)*0.25 UMS 7 

IFtKCM.EQ r . U)XP I' r CH-XI_F7 + (XTEZ-Xtm*0.'?0 IJMST 

I F ( KC M • EO •' r >IXPI rr H ^XLFZ+{ XTEZ-XLEZ ) *0. T’5 IJNST 

I F ( KCN «FQ • 16) X°T 7 TH=XLEZ+( X T EZ— XLEZ ) *0 * 350 UN *3 ^ 

!F(KCM.FD. l7)XPITCH=XLF/+( XTE7-XLF7 ) *0. 1875 UMST 

!F(KrM.FO.’B)XPITCH-XLEZ*+(XTEZ-XLEZ) *0.4500 UNS7 


UN I MAG= < 0 • ♦ l . V UNS 7 

C C M EG A= ( R R FD+UN I M A G*C P E D ) *UMACH/R ETA 


ALPS P C= •* L F m C *3.1 A 159/180. 

$ I MARCH'S I NT AL n RBC ) 

COS APC=COS ( AL FPBC ) 

HO 6 NN= 1 , N r 
A A ( NN ) ~0 • / 

INSERT THE KPCNECKER DELTA (SWITCH FOR THE LOWER SIDES OF DIAPHRAGMS ) 

DO 9 I S = l, NS 
NXYS=NXY( IS) ; 

DO 9 1 1= l » NXY S 
l = TT+NT( IS.) 

J = T "9x9:9 

LOWER STDFS OF DIAPHRAGMS 

I F ( ( KNC P M,L ( I S ) • EQ .- 1 ) . AND • ( I S.GT.NSBOQY) ) J= I -NX Y ( T S ) 

NNN-IX J-U«NTf5TAL 
A A f NNN ) = l • 

C ONT INtIF / ;■ : ■■■■■■■: 

CONST -- 7 ./•» . 1 AI 59 ' v'^i V ■: 

DO 10 1 = 1* N TOTAL r 

S U M DR ( I ) = 0. ■ : ‘V - ' ■’ 

SOURCE ( I ) =0. . ■ ~ 

S C T R AN ( I ) =0 . - : ; , 

siGNPTi i ) =i .o -T:-:'. : , ■. ■ . .■ 

DO 180 J ~ A C E = l » N S v "• 

J S B T Q T = N XY ( J F A C E ) ; • 

DO 1.80 JJ=l» JSBTCT ' . 

J=NT { J C ACE ) + J J 


DO 1002 K=l, 3 

>AAl(K,1. ) = 1.5 A ( YPP ( K * J ) — YMP ( K ♦ J ) ) 
A A 7 (K,p )= 0.5*(YPM(K,J )-YVM(K t J) j 
AA2(K*1)=0.5*(Y°P(K ,1 )- YEW ( K » J ) ) 
A A? (K * 2 ) =0.5* ( YV P ( K » J )-YMM(K, J )) 
AT IK,.* )=AAl(K,l ) 

A l ( K 1 2 ) =A A ll K ♦ 2 ) 

A 2 ( K » l ) = A A ? ( K * l ) 
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■y 

l 

f 


1002 

inn 3 

c 

’ 004 
C 


1006 

C 

100 7 

C 


c 


1008 

c 


c 


tot 0 

c 

c 


A2(K,2)=AA2(K,2 ) 

CONTINUE 
on 1003 1-1*2 

AU 11 L)=DCIP^0( A l ( 1,1.) , A 1 ( Z,L ) , A l ( 3, L ) , Alt l ♦ L ) , A l ( 2 , LI , A l ( 3 , 1 l ) 

A 2 A2 ( L ) = 001 PRG( A2( l ,0 *42(2,0 , 42 ( 3 , L ) ♦ A 2 ( l ,1) , A? (2,1) , A2 (3,1 ) ) 

00 1004 K= 1 , 3 

IF ( A1.AK U.£0.0.)Al(K,l)=Al(K,2) 

IF( A1A1(2) .B0.0.IA! (K,2 )-AMK,l) 

IF(A2*Z( 1) .EQ.0.)A2(K, i)=A2(K,2* 

!F(A2A2(2).E0.0,)A2(K,2)=A2(K,1) 

C O NT 1 N1J F v 

no 1.006 1 = 1 , 2 - 

Al.Al(U=OCTPPr( Air ,1 ) , A’ (2 ,1 ) , At (3, L ) , at (1,0, A’ (2,1) , A’ (3 , if) 
A2A2(L) = DCTPR0( A2( 1,L),A2(2,L) ,?2(3,L) , '2(1,0 ,42(2,0 ,42(3,0) 

no 1007 l= i , 2 - - '£■ 

no 1007 M= 1 , 2 

A 1 4 2 ( L , M ) =001 05 c ( A 1 ( 1 , L ) , A l ( 2 ♦ L) , A 1 (2 , L ) » A 2 ( l , : M ) , A 2 ( 2 , M ) , A2(3,M ) ) 

41 C «? A 2(1 ) = S Q R T ( At AT ( ? ) * A 2 A 2 ( L ) - A 1 A2 ( 2 , l ) ** 2 ) 

•410542 ( 2 ) =80° T ( 41 4 1 (1 1*4 2 42 (1 ) -A 142(1, 1 )**Z) 
AlCRA2(3)=S0 D T(ATAl(t)*A?A2(2)-AlA2( 1,2) **2) 

A1CR72( 4)=SGR t ( A?. A I, ( 2 )* A2 42 1 2 ) -A1 A2 ( 2 ♦ 2 ) **2 ) 

on loop k = i , 3 i; ■ I 

AVAUKH0.5*{ A’ (K, T )4-A? (K,2) ) ' • tl| 

A V A 2 ( K ) =0 . 6* ( A 2 ( K , 1 ) 4- A 2 ( K ,2 ) ) 

Y NORM (1 ) = AV AT 12 )*AV4? { 3 ) — A V A ’ ( « )* A V 42 ( 3 ) 

Y NORM (2 ) = AV A 1 (31* AVAR ( 1 )-4VA 1 ( i)*'V A 2(3) 

YMoom ( 3 ) = A y a 1 U )*AVA2 ( 2 )- AVA 1 ( 2 )*A VA2( 1 ) 

SN( l)=YNnPM{ 1)/RFTA 
SN( 2)=YN0RM( 2 ) 

S M C 3 ) =YNQRM ( 3) 

ASN = SORT ( S,N( * )**2+SN( 2 ) **2+54' ( 3 )**?. ) 

A YN = SQPT ( Y^'ORM ( 1 ) «3«2+YN0RM( 2) ^'32+YNCRN ( 3) ^*2 ) 
no mio k= 1,3 

U N ( K ) = Y N 0 R M ( K > /AYN 
SNIJNC K)=SN(K ) /ASM 
CONTINUE 


SUNnNO=SUPPRn(UN( l ) ,IJN( ?) ,UN( 3) ,<JN( l) ,UN( 2) ,UN(3) ) 


no 17 0 1 F ? r E= l , K R 

I F ( K D I A F ( I FA C E ) . N E . KD I AF f JF AC E ) . AND . NSROD Y . LT . NS ) 0C TO 170 

I S R TO T= NX Y ( IF AC E ) 

DO 160 I I = 1 , 1 SR TOT 
I -NT ( I FACE ) ♦ I I 


ORIGINAII 
OF POOR QUAIITB 


v'V'l\v 

.1 •'.-l- h - 



00 160 ISYMMY-l ,NISYMMY 
DO 160 ISYMM7-l,NSY^MZ ; 

SI R MPT ( 2 ) = 3 * l SY Y my 

SIGNPTI 3) =3.-2*ISYMMZ . 

00 : 1101 ICQPN°=l»4 ' 

DO 1101 ICCOP-=l-».? 

SIGNVK TCCR.NR » I COOP )=0. 

SIGNV2( ICQRNR, ICOOR J=0. 

00 1102 K= l » 3 

PZ f K ) =PC ( K t J ) - P C ( K » I ) *SI GMPT ( K ) 

P700P Z=DOT PR 0 ( P 7 ( " lt°?( 2)»PZ(3) » P 7 ( 1 ) , PZ ( 2 ) , P Z I 3 ) ) 
QDOTUN=DOTPRG(UN( 1) ,UN< 2) » UN ( 3 ) » D Z( 1 > , P 7 ( ? ) , PZ ( 3 ) ) 
IF ( I.EQ.JJGO TO HOT 

I F ( (ABS(OCnTUN) /SORT( PZDOPZ ) ) . IT. 1 . OE-4 ) qOQTUN-O. 
CONTINUE 

on ’7 to K=T 

0 ( K t l ) = Y P M ( K * J ) - nr { K t I ) * F I G N 0 T ( K ) 

Q ( K t 2 ) = Y°P < K , J ) -PC ( K » D * S I GNPT ( K ) 
0(K,3)=YMP(K,J)-PC(K, n*SIGNPT(K) 

0 ( K , 4 ) = Y MM ( K , J ) -PC { K t I) * S I GMP T ( K) 
continue ■•>■■■ 

00 nil K= 1*3 

QM ( K , 1)*0.5*'(G ( K, 1 ) fO ( K ♦ 2 ) ) 
0M(K»2)=0.5*(C(K,2)+O(K,3) ) 

QM ( K » 3 ) =0 . 5*( Q ( K » 3 ) +Q ( K * 4 ) ) 
0M(K,4)=0.5*(Q(K,4)+0(K» 1) ) 

continue 

DO : V’P 2 K-!:,3 

K°I-K+i 

KP2-K+2 

1 F ( KPf ,.G T .^ )KPI=KP*-?. 

I FIX p 2 . G T . 3) K P 2=K P 2- 3 


QMCP A 2 ( K , 2 ) =0 M ( KP l r 3 ) * 4 2 ( KP 2 » 2 ) -QM ( KP 2 » 3 ) *A 2 1 K P l , 2 ) 
CONTINUE 

00 1120' ICCPNP-If 4 . 

DO 1 12 D !C0GR = l »2 
M DEBUG (TCPPNP , ICOCR )= 0 

MK00E=0 

00 ?.!.30j I r nPN° = ’»4 - . • . 

KQDE ( ICORNP ) =0 
KCfilNT ( I CO°NR ) = 0 
KOOESt ( ICCS.NP ) = 0 
K0DEF2 ( I CORN 0 ) = 0 
K00ES3I I CORN 0 )=0 

OSIJPOI If OR NR ) - SUPPP C( 0 ( 1 , If CP NR ), 0 { 2 ♦ I CCRNR ) , 0 (3, I COP NR ) 


objsbs^I^ 


I Ofl, ICCRNP) ,Q(>, ICCRNR) ,0<3, ICnRNPn 

I F I OS UP Q ( I C 0 P N D ) .LE .0 .TP . G 1 1 , 1 CCP NR ) . GE.O >KOl)E l ICORNP 1=1 
MKQDF=MKOD£+KODE( T TCP NR) 

U30 CONTINUE ...i : - 

C ■ ' : ; 

TF(MKODE.ED.O JGO T 0 2209 i i ... . 

C ; ' 

C FOLLOWING IS TP CHECK THE DOUBLE INTERSECTIONS ALONG FT & DIRECTION A^O 

C TO EVALUATE T (-E CRITICAL VECTORS OF Q , AI «N0 S!GNVKI,2) A T MY C QR N 

0 POINT OF WHICH TEE ELEMENT IS PARTIALLY OUTSIDE THE MACH FOR ECONE 

C 

NK0DE=0 

r - : • -I ; ' : . 

nq 7220 I A 3=1.., 2 V; : ; : , 

I S I 0E= ( I A?- 1 »*?+? 

IFIKODE! I c I DE ) »FQ , 0 .4 NO . KODE ( l S I DE + 1 ) . EO . 0 ) GO TO 2220 
no 2210 K= l » 3 
22? 0 A A A 2 ( K I = A A 2 ( K , I A 2 I 
C 

CALL CP VEC T 1 1 SIDE » AAA 2, A AT , 2) 

C 

2220 CONTINUE 

r 

C FOLLOWING IS TO CHECK THF OOUBLE INTERSECTIONS ALONG CSl DIRECT I ON AND 

C TO EVALUATE THE CRITICAL VFCTCPS OF 0 , A 2 AND SIGNVK 1,1) 

r 

KODE! 5 ) = KODE ( l ) 

DO 2.24.0 141-3,2 
TSIOE=( IU-1JA2+2 

•IF (KODE! ISIDE) .FO.O.ANn.KOOE! I S I DE + l ) . EQ. 0) GO TO 2243 
DO 2230 K=T,3 
2230 AAUIKIsi.ntK, I.M ) 
r 

C AL L CR V EC T( ISIDE, A A At , A A 2 , 1 ) 

C 

2240 CONTINUE 
?299 CONTI NUF 
C 

I F ! MKODE »Ep . 4»AND »NKODE» EC.C ) GO T 0 160 | 

C . j 

IFtKODEC KEC.O.QP.KnDEm.FQ.C.OR.KCBINTU) .EQ.DGO TO 3310 | 

KODES 2 1 1 ) ■= 1 | 

KODES 2 ( 2 ) = l ■ ■' ; .'...I...;: \ j 

33LO I F ( KODE ( ; ) . ED . 0 « OR . KODE ( 4 ) . EQ . 0 « OP. . KOB I NT ( 3 ) . EO • 1 1^0 TO 3320 
K0PES2 ( 3 ) =1 
KODES 2! 4 ) = \ 

332) ! F ( K n DE ( , EO * 3 . OP « KODE ( 4 ) . EO « 0 • OR . KCB I NT ( A- ) .ED .MS 3 T 0 33^0 

KODE si ( I ) = 1 

KOOES l ( 4 )= 1 • 

333) I F ( KPDE ( 2 ) . EO . 3.CP .KODE t 3 ) . EQ • 0 .OR • KDR TNT ( ? ) . ED .1 ) GO TO 3340 j 

KODE c 1 ( 2 ) =1 
KODESl!? ) = 1. 

3343 I - ( KODE! I ) .EQ . I . AND .KCDE ( 2 ) . E 0.3) K ODE S3 It;) =* 

IF! KODE ( 2 ) . EO . 1 . AND .KCDE ( i 1 . EC *0 ) KODE S3 ( 2 ) = l 








!F(KPDEt 3 ).E 0 *T.AND.K 00 E( 4 ) .EQ.O)KQOES 3 ( 3 ) = l 

I F( KGOF( 4 ) .EC, l . ANTUKODEC 3 ) .EQ.O) KQ 0 ES 3 ( 4 )=l : 5 : . 

I F ( KHB I NT ( u .NE.UGO TO 3350 
• KODES 3 CT )=' 

KQ 06 S 3 ( ? ) = l 

3350 IFf K 0 BINT( 3 ) .NE.l IGO TO 3360 
KHDES ' 1 ( 3 ) 

KODES 3 ( 4 ) =1 :v-'~ 

3360 CONTINUE 
C . 

c ’ c=0 * ‘ 1 JNST 

° Z S 0 P Z = S U P P R 0 ( P Z ( * )t D Z( 3 ),PZ( 3 )fPZn ),P 7 .( 3 )»PZ( 3 i) UNST 

TF(°ZSP p Z«GT . 0 . iPr=SO p T{ 0 ZSPPZ) UNST 

CnSHU=(CEXP( C€^EGA*PC )+CEXP(-COMPGA^?.C) )^ 0.5 
S I N H U = ( f F X P ( r C M E 0 A * o C ) - ** E X P ( - C 0 V F G A * P f. ) >* 0.5 
CMoxDe=rOSHU-CQ y FGA*P :*SINHU 
CMPXSC=COSHU*CEXP(COMEGA*U«ACH*PC.( i»J ) ) 

C 

00 155 I CO^NP -l ,4 

GO TP, ( 5 502,5 50 4, 55 06,5 503 ) , ICOPNR 
5532 CONTINUE 

S ?GN 12 ~ - 1 . ' 

ICSI=i C; 

!?TA =2 - -'-Tc = /£ Cy<,- V v -C? 

GO TO 5510 
5504 CONTINUE 

SIGN 12 = + 1 . . 

ICSI = X 
teta=i 
GO TO 5510 
5506 CONTINUE 

$ IGNT 2 =- 1 . 

■. ■ - : res 1=2 •. . - 

Go TO 5510 
.,5503 CONTINUE 
\ SIGN 12 =+ 1 . ■ ' 

I ICS 1=2 

. [ '■ TE T A =2 ■ ■: '■ 

5510 CONTINUE 

' : fr, -C-V.. ■- : v ■■■■ 

00 5520 K= l ♦ 3 

QV(K)=Q(K, iCORNR) ; • CC 

■ • aiv(k)-a*. ( k, i eta ) i ■■■ ■,/ ■"'•■■.c- 

. ■ a 2v (k )=a 2 ( k , ics i ) ’■ T CT 

OCR A 1 ( K ) =OMC~ A 1 . ( X » I ETA ) ■ 

' 1 QC c A ? (K )=QWOo tft (K , ICS I) * 

.. 5520 . CONTINUE. "■. :-yC C V C,;;/' ' 

: C 1 ,.c--v;.c ; 

QQ=SQRT( A B S ( Q S U p Q ( I C PR NR ) ) ) 

. C ::-CU C.;..; v- C -C C 

A LOG 1 = 0 . '■ i.'V T t- v - : *■ ■■■■C'V, 

1 F ( KOOES l ( ICORNR) .FQ. I ) GO TO 5530 TCvVy — ■ / 

’ ■ CALL LGOt ICORNR, A IV, OCR Al ,AL 0 G 1 , l ) 





5533 CONTINUE 
C 

AIG&3=0, .'f‘| ' V ; ' 

IE|KnOES2( ICTRNR) .EOVtlCf Tr 5540 
CALL L CG I 1C G R N ® » A 2\i5> Q C'R A 2 » A LO G 2 » 2 ) 

5540 CONTINUE 

C- ■ .... ;.]■ 7 ; . V..r . =•--' / V v/. - C 

TAN 3=0. 

I F ( ODO T UN • EO . 0 • ) GO rr ) 5550 
IF(KOC c ( ICCRMP ) ,p«?„ l) GO TC 5545 

HNUMER = - SU3 or 0 ( QCP AT ( T) f OCR AT I 2 ) » QCP A 1 ( 3 ) » OCR A2 ( 1 ) • OCR A 2 1 2 ) 

1 OCR A 2 (.3 1 ) 

0EN0F=CQ*QnCTUN^AlCRA2( I CORN? ) : • . 

I F( OENOM.NE .0,. ) T i N P - A T A N ° IHMJNER , OENCNJ 
GO TO 5550 

5549 CONTINUE 

T F ( KODES a ( ! C CR NR ) , EG • 0 ) GO to 5550 
IFMDEEUOI ICORMP ,21 .FO.OJF ALL CE n UG(505) 

SONAR G=S IGNVl ( ICOPMR , 2 ) *SIGNV2 ( ICQRN 0 , 2 )*0D0TUN 
I FI SON ARO.L T . 0 . ) S GNTA N = - * , 

I F ( 5 GNAR G .EO.O. >SGN t A N = 4-0 . 

IF ( SGNARG . GT . 0 . )SGNTAN=+ l • 

T A NP= SON T A N * ' .573795 

5550 CONTINUE ■ ; - Cy;-;. 

COEFFl=SUPPRQ(UN(i) tUNI 2) ,IJN( 3) »QC p Al M ) » QC P A 1 ( 2) t OCR A l ( 3 )\ ,G 
COEFF2=SUOPRO(UM( 1 ) »UN( 2) ,UN< 3 ) »QCP A2I 1 ) ♦ QC R A 2 1 2 ) » 0CRA2I 3 ) k, <) 

c ■ ■ ' 1 . 

SRC!NT=-CONST*S IGM 2/ SUNONO* 

1- (-COEFF1*ALOG1+COEFF2*ALOG2-QOOTUN*TANP) 

c 

D8T INT =-C0N$T*S ION’ 2~TANP 
C 

SON I NT= l . 

I F ( I S Y M M Y . E 0 » 2 1 S 0 N I N T = S G N I N T * K S Y M m y 
I F I I S YMMZ . E Q . 2 ) SON I N T =S GN I N T * K $ Y M M Z 
C 

I F I J F A OF. G T • N S 3 CD Y ) GO TC 7.48 
C 

C AIRCRAFT 

C 

NMN=I +-( J-T )*NTO T AL 

c . 

A A ( N N N ) = A A I N N N ) - $ G N I N T * D B T I N T * C M P X 0 B UN ST 

C . '■ 

ARG = A&5( (PCI 2 ♦ J I-BODYR ) /( C. 5*<:P A N- n 00YR ) ) 

IF! I ZZZ. BO. lM ) XNIJMST=-I '’REn+UNIMAG’S'CREO )*SNUN( 3) * 

1 I ).-,0 343*APPFC . 7 0 ’ '5 5 * A® G^* 3- 1 3 4 3 8 =«« A R £**5 + 0. 25*37* ARG^ 4 1; . U>1 3 

IF! ! 7 Z7 .GE. 20Q . AND • IZ ZZ .L T . 300 ) XNUN S T =+ I 00 FOfUN l MAO’S'GR ED ) *SNIJN( 31 
IF( IZZZ.GE.300.AN0. UNST 

IXNUNST=-K (RPEOfUN I ^AO^C RED) *(Prr » J ) *BETA-XP I T CH ) +C. .0 ) ♦SNIJN l ? ) 

C / FOLLOWING IS FOP. 

■/ I c ( I Z 2 7 . LT .500 ) GO TC 499' 7 y: ; 'y ; - ' il ■■ 





X NUMST=+ ( ( RR'e 0 >UN I m AG*C R E 0 ) * { PC ( t , J I * 3 FT A-X P I TC H ) ♦ 1 • 1 ) * SNUNI ? ) 
XNTR AM=+ ( o ? EQ4IJ\ I PA G^C^ED )* SNUNI'* ) ' ' 

SCTRA N( I )=SCTO AN { I ) +SGN!NT*$RC I NTT X N TR A NT CM PX S r 
CONTINUE 

s nuocem = source i 1 1 *$ gnint* srci nttxnun st*c m px sc 

GO TO 150 * 

’48 CONTINUE , V 

r v . . 

C DIAPHRAGM r j 

C ■■■:.).. : r ■■ : :■■■ 

TFIKSYMMZ .NE.OGQ tq 149 
JPHt= J 
JCH I = J 

IFIKDt AF ( T FACF ) .EO.-L > JPH 1-= J-NXY < JF A CE I 
IF(KDIAF( IFACE) .frO. + l ) JCH I = J + NXY ( JFACF ) 

N° H I = 1 + ( JP H I - 1 ) TNTCT A L 1 

N CM I = t + ( J C H I - I ) *N TO T A L 

c 

AAINPHI J=AA(NPHI )-SGNINTTD8TTNT*CMPXDB 
AAINCHI ) = AA( NCHI ) -SON I N T * S* C I N T *KD I AF ( I FACE > 
r 

GO TO 150 

149 TFIKSYMMZ.EO.-- . )GO T 0 1491 

r 

C SYMMETRIC : CHI=0 

c ' ■ . 

N PH I - I + ( J- 1 ) T NT C T A L 

A A ( N P H I )~ A A ( N P H I ) - S GN I N T T 0 P T t N T TC M P X 0 R 
GO TO 150 
1491 CONTINUE 

c ' ; ; 

C ANT I S VMM E TO IC : PH!=0 

C 

NCHl = I + ( J-DTNTOTAL 

AMNFHI ) -'\* (:NCHI )-SGMTNTTSRCTN t TCMPXSC 

150 CO NT I NU r 

IF! 1.01.0 )GO TO * 55 
C I F ( MKQDE .£0.9)00 T 0 155 

WRITE (6, 66) I , J , ISYMW, I SYP-MZ , ICORNP ,KODE ( I CORNS ) » S RC I NT , 03TINT 
1ALOGT ,rOFFFl,MCG2*CnEEFT 

66 F0PMATI1X.6I4, 7E15.6) . . 

155 CONTINUE 

16 0 CONTINUE . • > ■■■■ 

170 CONTI MU F V" 

130 CONTINUE 

return’ ' 

END ■ ■ ■ .. Uv : >. 
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SUBROUTINE CR VEC T ( l SIDE » V l » V2 » ! COOR ) 

C OMMDN/7 7 Z 5 6 / NKODE , KD BINT (4 ) t MDEBUG ( 4, ? ) , QM ( B f 4) , 

C CVMON/Z Z Z77/ QV| ? ) »CO ? KCDE ( 5 ) » S I O NV ", l 4*2 ) * SIGN V’ ( 4* 7 ) 

C CM MON/ 7 7 ZS3 / I » J ♦ I SV'IMY » I SYMMZ 

DIMENSION GCQQR ( 2 ) , C° V2( 3 ) *OQC (3 ) ♦ QQCXI 2 ) t Vi ( 3 1 * V2 ( 3* 7 ) 
SPVIV1=V1( 1)«*7-VK 2)**2-Vl(3)**2 
I F ( $ P V l V l . E Q . c . ) 0 C TO 1999 ■ 

S D 9ZQ/.=0V( 1, IP tqe )**2-QM( 2* I S IDE)**?-Qv( 3 , 13 1 OF ) **2 
S D V10?=VI ( 1 ) * CM { 1 i ! S I OF ) — V 1 1 2 ) *QM( 2 , 1 SIDE )-Vl I 3 ) *QM( 3 , ! S I DE ) 

^v;;.:.,/;^:&|SCP.M*SPV^pZ**2-S.P.az^Z^S'l>Vt■Vi■ 

I FIDI SCP w.LE.O. )GO Tq 199 c 
AvcnoR=-s D vnz/s°vivi 

X&VEP G=QM( ? f ISIDE )+AVCCQP*Vl III 

SQDIS=SQRTCDISCRM)MR S(S D V1VI) .V 

. GCOOR (1 ) ==AVCODR + SQD IB 

GCOOR ( 2 I = A V r OCR - SOD I S V ' 

oqgxi i )=qm( l» igioei+gcqod (nn’iiii 
OQCXI 2l=OM( l, I3IOE)+GCOnRI2)«Vl(l) 

IFI (QQCXI? )*QQ0X{?.) ).GE.Q.)GO TO 100 
: GCOQP ( 1 1 = 4 VC 00 P - SOD I S ■ * ■ 

GCOOR I 2 ) = A VCOD R + SQD IS 

QQG X ( t J=QV{’ , I S IDE ) +GCQCR ( 1 >*V? (I ) 

OQCX ( 2 ) =0M( l , IS IOE)+GrOOR (2 JT-V1 ( 1 ) 

100 CONTINUE 

DO 300 I ENO = l , 2 

I F( OCOOR I I END ) .GT . 1. . OR .GCOOR IT END). LT. -I. .OR. OQCX 1 I END) #C>T.O. ) 
1G0 T 0 300 

IF! IS IQF.EQ, ]. ) TC'~PNR= ISIDE +I2-IE NO) 

I FITS IDF. EQ. 2 ) I COR MR = I S T 0 E+I I END- I ) 

I FI TS IOE . FO. 3 ) I COR NR- I S ICE+I I END- l ) 

I F ( IS I OF . ,~Q . B ) I C OR N p = I IS I D E + 1 I END— 2 ) #3 
I F ( MDEBUG I TCC C NP V tCCOR ) . EC .1 ) CALL DEBUG I 5 10) 

I F ( KODE ( I CORNS ) . EQ.C ) GO TC 300 
DO' 200 K«l, 3 

CRV2 1 K ) -0 . 5# I V2 ( K * l) + V2(K,2) ) EO. B*0Cf)0P ( TEND ) * ( V2 I K 1 1 ) -V? I K » 7 ) ) 

' 200 QQO (K ) =0M { K 1 1 S I HE ) +GCOOR ( I FND )*V! (K ) 

QSUPV l-QOC ( l) *V II 1 ) -OQCI ? ) * V 1 ( 7 )—QQC { 3 ) *Vl( 3 ) 

QSU D V2=0QC( I )*CRV2( 1 ) -QCG I 7 ) V2 ( 2 ) -QOC ( 3 ) *C 7 V2 I 3 ) 

SIGNVK I r 0° N° » I r OCP )= CSIJP V l / AB S (QSlJ<° VI ) 

S IGNV2I TCQRNR , IOOOP ) = QSURV2/ABS(0SURV2 ) 

M DEBUG ( I CCPNR ♦ iqoOR )= 1 : 

300 CQNTIMUF . '■■■ ' J 

I F I KODE! I S ID^ ) . -Q .O.OR. KODF ( I SIDE* l ) . EQ . ) )G0 TO ) 999- 

IF( AVC00R.GT. I.. OR. AVCDOR.LT. -V.. OR. XAVERG.GE.O.. 'OR. SPOZOZ.Lr.O. ) 

1G0 TO 190 9 ■ 

• KDOINTI IS!DF) = 1 -v ... 

■ NKOOE=* ;■■■ 

1909 CONTINUE . V. v : d ! V ; 

RETURN ' . iV- • - ■ \/" ! " :■ V v- 

END .'.'.'.i ■ -J ■ 
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SUBROUT INE LOG ( ICCRNR » A V, OCR A V , ALOG » I CGCR ) 

C CMMOM/Z 7. ZR7/QV (?) , 00 rKODE( 5 ) ♦ SIGNV7 { A* ? ) , STGNV2I 4,2) 
DIMENSION *V(3 ) ,QGR A V(3) 

e: ; v , \ : : Y ■ • ■ 

SUPDPCHXl ,Yl, Z1. , X2,Y2 tZ2)=X**x?-Y , ,*Y2-Z’ *V- 

c : ■ I r y ■■ ■ 

AVSP^ V=SU°PP f j ( AVI 1 ) ,AV( ? ) , A V ( 3 ) , A V < 1 1 » A V ( 2 ) * A V ( 3 ) I 
AVMQPMsY. 

I F ( AV S P AV . NF . 0 . ) AVNE 17 S 0 R T ( A 0 S ( A V S ° A V ) ) 

I F t RODE ( T CORN P ) ,■ E C * I ) GO T Q 10 * 

QSU rJ AV=SUP rT Rn (0V( 1) ,0 V l ?) ,Q\n 3) * AVI 1) » AVI 2) , AV( 3 ) ) 

OCR ASP=SU P °* ? n ( OCR AV ( * ) ,QC P AV( 2) ,QCRAV<3) , OCR AVI L ) ,CCRAV< 2) , 

- 1 - j . qcp.av{3) ) 

IF1QCRASP.GF.0. )CALL DEBUG (520) 

0XA = S0PT(-3CP AS p ) L V 

I : F t A V'jSP ? V . GT . 0 . H IPG*> S T NH ( QQ* * VNO» Y/CX * ) /A-VNORS-. 

I F ( OSU p A V . L T. 0 . ) A Ll 1G=- ALOQ 
I Ft AVSP AV • EQ » C • ) A LCG~ CC/G SUP A V. 

I F ( AVSP AV . L T . 0 . ) 

IAL0G=-ARS IN (OSUPAV/SQRT ( QSUPA V*0SU p 4V-QQ*0Q*AVS°AV) l/AVNORM 
P. ETU C M 

10 C ON T T NUF 

I F I Avspav.gf.o. )ALGG-G* 

l = ( AVSP4V.I T . ) . ) A L 0G= - S I G N V 1 U C CPNP , I COPR) *} 57 379 5 / AVNFP M 

RETURN 

END 
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